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FOREWORD 


The  choice  of  functions  to  be  approximated  for  this  report  was  determined  by  the 
requirements  of  a  project  for  the  computation  of  surface  waves  in  the  wake  of  a  ship. 

The  computation  of  rational  functions  was  programmed  originally  for  the  Naval 
Ordnance  Research  Calculator.  The  development  of  rational  functions  was  pushed  to 
a  sixteen  decimal  digit  accuracy  before  the  destruction  of  NORC  would  make  all  this 
programming  useless  for  further  computation.  A  sixteen  decimal  digit  level  of  accuracy 
was  selected  because  this  is  the  level  for  double  precision  computation  on  the  IBM 
360  computers.  It  is  more  accurate  than  necessary  for  single  precision  computation 
on  the  CDC  6600  computers,  but  no  attempt  will  be  made  to  readjust  the  approximations 
back  to  a  lower  level  of  accuracy. 

Computation  of  constants  and  programming  of  checkouts  were  contributed  by  Mrs. 
E.  J.  Hershey  and  by  Mr.  W.  H.  Langdon.  The  manuscript  was  completed  by  11  July,  1971. 
The  report  was  reviewed  administratively  by  Mr.  J.  H.  Walker,  Jr. 


Approved  for  release: 

R.  I.  ROSSBACHER,  Head 


Warfare  Analysis  Department 
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ABSTRACT 


A  few  special  functions  have  been  approximated  in  the  complex  plane  with  rational 
functions.  The  error  bounds  of  the  approximations  conform  to  the  Chebyshev  criterion. 
The  rational  functions  have  been  converted  into  equivalent  expansions  in  terms  of  the 
singular  functions  for  poles  and  residues.  Analyses  and  programs  are  described  for 
ten  functions. 


INTRODUCTION 


The  choice  of  method  which  will  be  used  in  the  preparation  of  a  special-purpose 
subroutine  is  dictated  by  the  requirement  for  efficiency  in  the  specific  application  for 
which  the  subroutine  will  be  used.  The  choice  of  method  which  will  be  used  in  the 
preparation  of  a  general-purpose  subroutine  is  dictated  by  the  requirement  for 
versatility  in  a  computer  system.  A  sensible  policy  is  to  prepare  a  subroutine  in  such 
a  way  that  the  rounding  error  in  output  is  compatible  with  the  inherent  error  in  input 
for  any  single  precision  computation  within  the  capacity  of  the  computer.  With  such 
a  policy  it  is  likely  that  the  subroutine  will  have  the  greatest  range  of  application. 
The  effort  which  has  been  expended  in  its  preparation  will  have  the  greatest  reward. 

A  number  of  special  functions  can  be  expressed  as  the  sum  of  an  ascending  power 
series  in  the  argument  or  as  the  sum  of  a  descending  power  series  in  the  argument. 

The  ascending  power  series  is  absolutely  convergent.  The  truncation  error  in  a  finite 
sum  can  be  made  less  than  any  arbitrary  bound  by  taking  enough  terms  in  the  series. 
When  the  argument  is  outside  a  circle  of  unit  radius  the  magnitudes  of  the  terms 
may  increase  at  first  with  increasing  order  to  a  maximum  magnitude  although  they 
always  decrease  ultimately  with  increasing  order.  When  each  term  is  computed  to  a 
finite  number  of  digits  the  rounding  error  in  the  sum  is  determined  principally  by 
the  rounding  error  in  the  largest  term.  In  that  part  of  the  complex  plane  where  the 
terms  all  have  the  same  sign  the  rounding  error  is  relatively  ineffective  because  the 
sum  is  larger  than  the  largest  term,  but  in  that  part  of  the  complex  plane  where  the 
terms  have  alternating  signs,  the  rounding  error  can  be  completely  fatal  because  the 
sum  is  very  much  smaller  than  the  largest  term.  The  ascending  power  series  is 
inadequate  outside  this  contour  of  limited  extent  because  of  rounding  error. 

The  descending  power  series  is  asymptotic.  The  magnitudes  of  the  terms  decrease 
or  increase  at  first  with  increasing  order,  but  they  always  decrease  to  a  minimum, 
and  finally  increase  with  increasing  order.  Where  the  terms  alternate  in  sign,  the 
truncation  error  in  a  finite  sum  is  less  than  the  smallest  term  included.  The  truncation 
error  is  less  than  any  arbitrary  bound  when  the  argument  is  outside  a  circle  of 
sufficiently  large  radius.  The  descending  power  series  is  inadequate  inside  this  circle 
of  large  radius. 

There  is  a  zone  between  the  circle  of  unit  radius  and  the  circle  of  large  radius 
which  cannot  be  reached  by  either  of  the  classical  series.  The  conventional  approach 
to  this  region  has  been  to  seek  a  converging  factor  which  can  be  applied  either  to 
the  leading  term  or  to  the  smallest  term  of  the  asymptotic  series  and  thereby  achieve 
a  satisfactory  level  of  accuracy.  Various  schemes  have  been  proposed  in  the  literature 
for  evaluating  the  converging  factor. 

In  the  method  of  Airey1,  the  smallest  term  is  factored  out  of  the  following  divergent 
terms  to  obtain  a  series  of  divergent  factors,  each  of  which  is  expanded  then  in  inverse 
powers  of  the  order.  Successive  terms  with  the  same  inverse  power  of  the  order  are 
replaced  by  their  Euler  sums.  The  terms  in  the  summation  converge  for  a  limited 
range  of  the  argument,  but  the  results  of  summation  are  rational  functions  which 
can  be  evaluated  at  any  argument. 

Although  the  substitution  of  an  Euler  sum  for  a  divergent  series  is  dubious,  the 
Airey  formulae  for  the  exponential  integral  and  the  probability  integral  have  been 
confirmed  rigorously  by  Murnaghan  and  Wrench5,6  through  transformations  of  variable 
and  integrations  by  parts.  They  have  extended  the  Airey  converging  factors  to  many 
more  terms.  Although  their  formulae  give  a  great  improvement  in  accuracy  when  the 
smallest  asymptotic  term  is  small  anyway,  they  are  not  especially  helpful  when  the 
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argument  is  near  the  unit  circle.  It  is  not  clear  that  the  terms  in  the  converging 
factor  are  not  themselves  asymptotic.  Cf.  Page  22  of  Reference  5. 

In  the  method  of  Stieltjes2,  the  integral  representations  of  the  Bessel  functions  are 
converted  into  double  integrals  whose  integrands  contain  a  binomial  in  the  denominator. 
The  reciprocal  of  the  binomial  is  expanded  in  power  series  with  remainder.  Integration 
of  the  power  series  term  by  term  recovers  the  terms  of  the  asymptotic  series  while 
the  remainder  is  retained  as  an  integral.  A  transformation  of  the  variables  of  integration 
and  an  expansion  of  the  integrand  in  a  power  series  provides  an  evaluation  of  the 
remainder  in  inverse  powers  of  the  order.  This  method  is  on  firm  ground  insofar  as 
its  errors  are  limited  to  the  approximation  of  an  otherwise  exact  formula  for  the 
remainder,  but  it  requires  the  evaluation  of  an  exponential  function  in  the  evaluation 
of  the  remainder. 

The  Stieltjes  method  has  been  extended  to  more  functions  and  to  higher  orders  by 
J.  T.  Moore3,4  in  this  laboratory.  Insofar  as  the  Airey  method  and  the  Stieltjes  method 
both  approximate  the  remainder  in  a  series  of  descending  powers  of  the  order,  it 
would  be  expected  that  the  methods  would  be  equivalent,  and  this  has  been  confirmed 
by  comparisons  in  Reference  27. 

In  the  method  of  Dingle9,10,  basic  functions  are  defined  by  the  Borel  sums  of  their 
asymptotic  expansions.  The  remainders  of  the  asymptotic  expansions  provide  a  set  of 
basic  converging  factors.  In  the  asymptotic  expansion  of  the  confluent  hypergeometric 
function  the  remainder  is  an  integral.  The  converging  factor  of  the  hypergeometric 
function  is  expressed  as  an  integral  of  the  basic  converging  factors.  The  integrand  of 
the  converging  factor  is  expanded  as  a  series  in  ascending  powers  of  the  argument 
and  is  integrated  term  by  term.  The  method  is  justified  insofar  as  it  reconstructs  the 
integral  representation  of  the  confluent  hypergeometric  function.  Dingle’s  method  has 
been  determined,  at  least  in  the  case  of  Bessel  functions,  to  be  itself  asymptotic  with 
significant  truncation  and  rounding  errors.  Cf.  Reference  27. 

In  a  study  by  Tonelli11,  consideration  has  been  given  to  the  uniqueness  and 
convergence  of  polynomial  approximations  which  conform  to  the  Chebyshev  criterion 
on  a  closed  contour  in  the  complex  plane. 

In  the  method  of  Hastings12,  the  converging  factor  would  be  approximated  by  a 
rational  function  with  an  error  which  conforms  to  the  Chebyshev  criterion.  The 
maximum  absolute  error  in  each  loop  of  the  error  curve  would  be  the  same  in  every 
loop. 

The  absolute  value  of  any  analytic  function  within  a  closed  contour  is  everywhere 
less  than  its  absolute  value  on  the  closed  contour.  The  error  is  analytic  if  the  function 
and  its  approximation  both  are  analytic.  The  rational  approximation  is  analytic  if  the 
roots  of  its  denominator  lie  outside  of  the  contour.  The  error  then  is  bounded  everywhere 
within  the  contour  if  it  is  bounded  on  the  contour.  If  the  contour  lies  within  the  region 
which  can  be  reached  by  one  of  the  classical  series,  but  encloses  the  critical  region 
which  cannot  be  reached  by  either  of  the  classical  series,  then  the  approximation  will 
span  the  critical  region. 

The  rational  approximation  can  be  specified  at  a  series  of  points  along  the  contour13. 
The  points  of  specification  can  be  situated  at  antinodal  points  or  at  nodal  points  in 
the  variation  of  error  along  the  contour.  The  antinodal  specification  would  give  a  lower 
error  bound,  but  it  would  be  difficult  to  control  the  number  and  location  of  the 
antinodal  points.  In  the  nodal  specification  there  is  a  maximum  of  error  between  each 
consecutive  pair  of  nodes.  The  nodal  points  can  be  adjusted  so  as  to  equalize  the 
maxima  of  error  between  nodes.  The  algorithm  which  was  used  in  the  development  of 
rational  approximations  is  given  in  Appendix  A. 
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An  effective  contour  of  approximation  lies  along  a  half  circle  of  unit  radius,  along 
the  imaginary  axis,  and  along  a  half  circle  of  infinite  radius  to  form  a  closed  contour 
which  contains  the  critical  region  of  approximation.  Although  the  error  of  approximation 
is  bounded  only  within  the  contour  of  approximation,  the  rational  approximation  is 
good  enough  to  be  used  in  an  extensive  region  outside  of  the  contour  of  approximation. 
The  rational  approximation  and  the  convergent  series  together  can  span  the  entire 
complex  plane  with  little  more  error  than  the  inherent  error  which  arises  from 
uncertainty  in  the  argument  itself.  The  boundary  line  between  the  regions  of  application 
of  the  different  formulae  is  determined  empirically. 

The  preparation  of  rational  approximations  with  single  precision  requires  source 
data  with  higher  precision  to  serve  as  foundation  and  standard.  The  achievement  of 
better  than  single  precision  accuracy  was  not  possible  in  the  case  of  Bessel  functions 
without  recourse  to  converging  factors.  The  achievement  of  better  than  single  precision 
accuracy  was  possible  in  the  case  of  other  functions  through  sixteen— point  Gaussian 
integration  with  double  precision.  The  numerical  constants  for  the  double  precision 
computations  are  available  in  the  literature16-20. 

The  rational  approximations  can  be  expanded  into  the  sums  of  singular  functions 
which  are  useful  for  further  analysis.  The  singular  functions  represent  poles  and 
residues  at  the  roots  of  the  denominators  of  the  rational  approximations.  An  algorithm 
for  the  conversion  of  a  rational  function  into  a  polar  expansion  is  given  in  Appendix 
A.  The  results  of  approximation  are  given  in  Appendix  B. 

The  errors  of  approximation  are  located  on  spiral  curves  in  the  complex  plane. 
Representative  curves  are  illustrated  in  References  25  and  26.  Maximum  errors  in  the 
polar  approximations  along  the  imaginary  axis  and  along  the  arc  of  unit  radius  are 
tabulated  in  Appendix  C. 

Justification  for  the  success  of  the  polar  expansion  is  to  be  found  in  a  paper  by 
Gautschi14  on  the  error  function.  The  asymptotic  expansion  of  the  error  function  is 
identified  with  the  asymptotic  expansion  of  a  Stieltjes  integral,  and  approximation  of 
the  Stieltjes  integral  by  Gauss-Hermite  quadrature  leads  to  an  expansion  in  poles  and 
residues. 

In  the  method  of  Goldstein  and  Thaler15,  the  Bessel  functions  are  evaluated  by  a 
cycling  of  recurrence  equations  from  such  a  large  order  that  the  Bessel  function  is 
negligible  in  comparison  with  the  Weber  function.  A  large  number  of  cycles  of  recurrence 
is  required  when  the  argument  is  large. 


CUTS 

For  a  general  purpose  subroutine  which  generates  a  complex  function  of  a  complex 
argument  it  is  desirable  to  restrict  the  range  of  argument  by  suitable  cuts  in  the 
complex  plane.  Otherwise  the  function  to  be  generated  would  be  multiple  valued  as 
the  argument  encircles  a  singularity.  The  call  lines  of  the  subroutine  would  be 
encumbered  with  parameters  which  would  specify  the  number  of  times  that  the 
argument  has  encircled  each  singularity.  In  the  present  system,  it  is  assumed  that 
the  complex  plane  is  cut  by  straight  lines  which  extend  from  each  singularity  to 
infinite  distance  in  a  direction  radially  outward  from  the  origin.  Whenever  the  argument 
happens  to  fall  exactly  on  a  cut,  it  is  assumed  to  belong  to  that  side  from  which  it 
would  come  during  a  counterclockwise  displacement  around  the  singularity.  When  the 
function  to  be  computed  is  for  a  path  which  crosses  the  cut,  then  the  path  can  be 
varied  to  encircle  the  singularity  and  the  function  can  be  incremented  by  its  change 
of  value  during  encirclement. 
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COMPLEX  GAMMA  FUNCTION 


Analysis 

The  gamma  function  r(z)  is  defined  by  the  equation 


z 

1  +  - 
n  . 

where  y  is  Euler’s  constant.  The  gamma  function  has  poles  at  the  negative  integers 
such  that  the  residue  of  the  nth  pole  is  (— l)n/n!. 

For  a  small  argument  the  reciprocal  of  the  gamma  function  is  given  by  the  Bourguet 
convergent  series  and  for  a  large  argument  the  logarithm  of  the  gamma  function  is 
given  by  the  Stirling  asymptotic  series.  Intermediate  regions  can  be  spanned  by 
recurrence  relations.  A  rational  approximation  is  not  necessary. 

The  gamma  function  of  an  argument  with  a  negative  real  part  is  expressed  in  terms 
of  the  gamma  function  of  an  argument  with  a  positive  real  part  by  the  reciprocal 
equation 


r(z)  =  —  n 


r(z)r(i  -  2)  = 


7T 

sin  7 tz 


(2) 


It  is  necessary  to  evaluate  series  expansions  only  for  arguments  with  positive  real 
parts. 

If  the  argument  x  +  iy  satisfies  the  inequality 

x*  +  yzH  1  (3) 


then  the  gamma  function  is  derived  from  an  ascending  power  series.  The  reciprocal 
of  the  gamma  function  is  given  by  the  equation 


1 

r(i  +  2) 


00 


=  E  cmz 


(4) 


for  which  the  coefficients  cm  are  listed  on  page  256  of  Reference  21. 
If  the  argument  x  +  iy  satisfies  the  inequality 

xz  +  yz>  32 


(5) 


then  the  gamma  function  is  computed  from  a  descending  power  series.  From  the 
equations  on  page  252  of  Reference  23,  the  logarithm  of  the  gamma  function  is  given 
by  the  equation 


n-  1 


log  T(z)  =  (2  -  |)  log  z  -  z  +  4  log  (2tt)  +  E 


(-1 


2m(2m  —  l)z2 


where  the  Bernoulli  numbers  Bm  are  defined  by  the  equation 


(6) 


4  2  f  1 

(2m)!  (22m  -  l)7T2m  to  (zk  +  l)2m 


(?) 


Summation  of  the  series  is  continued  until  there  is  no  change  in  sum  or  until  m  =  18. 
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If  the  argument  x  +  iy  satisfies  the  inequality 

1  <  xz  +  yz  <  32 


(8) 


then  the  gamma  function  is  computed  with  the  aid  of  the  difference  equation 

r(i  +  z)  =  zr(z)  (9) 

If  n  is  the  integer  which  is  nearest  in  value  to  x  and  if  n  satisfies  the  inequality 

\z  —  n|2  ^  1  (10) 

then  the  gamma  function  is  given  by  the  equation 

r(z)  =  (z  -  1  )--(z  —  n  +  l)F(z  —  77+1)  (11) 

for  which  T(z  —  n  +  1)  is  evaluated  from  the  convergent  series.  If  n  is  the  smallest 
integer  which  satisfies  the  inequality 

\z  +  n\z>32  (12) 


then  the  gamma  function  is  given  by  the  equation 

r(z  +  n) 


r(z)  = 


(13) 


Z  -(z  +71-1) 

for  which  T(z  +  n)  is  evaluated  from  the  asymptotic  series. 
Programming 

SUBROUTINE  CGAMMA  (MO,  AZ,  FG) 

FORTRAN  SUBROUTINE  FOR  COMPLEX  GAMMA  FUNCTION 


The  mode  of  operation  is  given  in  MO.  The  real  and  imaginary  parts  of  the  argument 
z  are  given  in  array  AZ.  The  complex  gamma  function  is  computed  by  series  expansions 
and  recurrence  relations.  If  MO  =  0,  the  real  and  imaginary  parts  of  the  function  T(z) 
are  stored  in  array  FG.  If  MO  =  1,  the  real  and  imaginary  parts  of  the  function  log  F(z) 
are  stored  in  array  FG. 
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COMPLEX  LOGARITHMIC  INTEGRAL 


Analysis 

The  logarithmic  integral  L(  1  +  2)  is  defined  by  the  equation 

L(  1  +  2)=  f  -log(l  +  t)dt  (14) 

Jo  ^ 

The  integrand  is  analytic  at  t  =  0  but  has  a  singularity  at  t  =  —  1.  The  complex  plane 
is  cut  from  the  singularity  to  along  the  negative  real  axis. 

The  entire  complex  plane  is  covered  by  a  pair  of  series  which  are  centered  at  the 
origin.  The  rates  of  convergence  of  the  series  are  slow  in  an  annular  zone  which 
contains  their  common  circle  of  convergence.  The  annular  zone  is  covered  by  a  pair 
of  series  which  are  centered  at  the  singularity.  A  rational  approximation  is  not 
necessary. 

If  the  argument  x  +  iy  satisfies  the  inequality 

xz  +  yz<\  (15) 


then  a  Taylor  series  expansion  of  the  logarithm  in  the  integrand  leads  to  the  series 


L(l+z)  =  -  £ 

m= 0 


(~2)m+1 

(m+  l)2 


which  is 
If  the 


convergent  for  2  inside  the  circle  \z\  =  1. 
argument  x  +  iy  satisfies  the  inequality 


xz  +  yz>  9 


(16) 


(17) 


then  the  logarithm  in  the  integrand  is  replaced  in  accordance  with  the  equation 


log  (1  +  t)  =  log  t  +  log  ( 1  + 
Expansion  in  powers  of  1/t  leads  to  the  equation 

7T2  00  1 

L(1  +  2)  =  -  +  £ 

771=0 


6  _~0  (m+  l)2(-2)m+ 


I  +  |log2Z 


in  which  the  constant  of  integration  is  derived  from  the  equation 

1  7TZ 

—  log  (1  +  t)  dt  =  — 
t  12 

and  the  series  is  convergent  for  2  outside  the  circle  |z|  =  1. 

If  the  argument  x  +  iy  satisfies  the  inequality 

(1  +  x)z  +  yz^i 

then  integration  by  parts  gives  the  equation 

J  j  log  (1  +  t)  dt  =  log  (-2)  log  (1  +  z) 


r  log 

Jo  1 


log  (-t) 


+  t 


dt 


(18) 

(19) 

(20) 

(21) 

(22) 
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Expansion  in  powers  of  1  +  t  leads  to  the  equation 


7T2  “  ('l+z')m+1 

Z(1  +  z)  =  -  —  +  Y,  +  log  (-z)  log  (1  +  z)  (23) 

6  m=0  (m+  l)2 

in  which  the  constant  of  integration  is  derived  from  the  equation 

P-1  1  tt^ 

J  -  log  (1  +  t)  dt  =  -  —  (24) 

and  the  series  is  convergent  for  z  inside  a  circle  with  center  at  —1  and  with  unit 
radius. 

If  the  argument  x  +  iy  satisfies  none  of  the  inequalities  above,  then  let  z  be  replaced 
by  u  in  accordance  with  the  equation 

u  =  —  log  (1  +  z)  (25) 

whence  the  integral  is  transformed  in  accordance  with  the  equation 


log  (1  +  f)  dt  =  — 


(26) 


From  formula  764,  page  90  of  Peirce’s  Table  of  Integrals22  the  integral  is  given  by  the 
equation 


u-  -  (— 1)"B„ 


i(1  +  ,)-.„+_  + £(sn+i)| 


u 


(27) 


where  the  Bn  are  Bernoulli  numbers.  The  series  converges  when  z  lies  within  a  pair 
of  contours  which  are  concentric  with  the  singularity  at  -1.  The  contours  are  defined 
by  the  equation  |ii|  =  2it  and  enclose  all  of  the  other  circles  of  convergence. 


Programming 

SUBROUTINE  CLGMCI  (AZ,  FL) 

FORTRAN  SUBROUTINE  FOR  COMPLEX  LOGARITHMIC  INTEGRAL 

The  real  and  imaginary  parts  of  the  argument  z  are  given  in  array  AZ.  The  complex 
logarithmic  integral  is  computed  by  Taylor  and  Bernoulli  expansions.  The  real  and 
imaginary  parts  of  the  function  A(1  +  z)  are  stored  in  array  FL. 
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COMPLEX  EXPONENTIAL  INTEGRAL 


Analysis 

The  complex  exponential  integral  Ei(z)  can  be  defined  in  terms  of  its  complex 
argument  z  by  the  equation 


(28) 


where  the  path  of  integration  lies  in  that  part  of  the  complex  plane  from  which  the 
positive  real  axis  is  excluded.  The  conventional  exponential  integral  of  real  argument 
x  is  the  real  part  of  the  complex  integral  when  the  complex  argument  z  is  real. 
Substitution  of  —t  for  t  gives  the  equation 


Ei(x 


Joo 

-x 


~Y~  dt  =  3ieEi( x  ±  0 i) 


The  sine  integral  Si(x)  is  defined  by  the  equation 

sin  t 


Si(x) 


sin 

Jo  ^ 


dt 


and  the  cosine  integral  Ci(x)  is  defined  by  the  equation 

cos  t 


Ci(x) 


poo 

J|*l 


t 


dt 


(29) 


(30) 


(31) 


The  conventional  integrals  are  recovered  from  the  complex  integral  by  the 
transformations 


t  -*  it 

and  are  expressed  by  the  equation 


z  ->  ix 


(32) 


Ci(x )  +  iSi(x)  =  ±  +  Ei(ix)  (33) 

where  the  sign  of  ±  if  is  the  same  as  the  sign  of  x. 

The  Stieltjes  form  of  the  exponential  integral -is  given  by  the  equation 

e~u 

Ei(z)  =  ez  - du  (34) 

Jo  z~u 

which  may  be  derived  from  the  substitution  t  =  z  —  u. 

If  the  argument  x  +  iy  satisfies  the  inequality 

xz  +  yz^l  (35) 

or  both  of  the  inequalities 

x2  +  yz  <  1600  -  x  +  0.064  y2  S  0  (36) 

then  the  exponential  integral  is  computed  from  the  ascending  power  series.  Substitution 
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of  —t  for  t  in  the  complex  integral  leads  to  the  equation 


Ei(z)  = 


1  —  e" 


dt  ■ 


It -j: 


df 

T 


1  —  e" 


df 


Evaluation  of  integrals  and  expansion  in  series  leads  to  the  classical  equation 

00 

Ei{z)  =  7  +  log  (-0)  +  XI 


v=i  m  m' 


(37) 


(38) 


for  which  the  constant  of  integration  is  Euler’s  constant  7  and  is  defined  by  the 
equation 

1 


7 


i 

Jo 


1  1  -  e  1  -  e  1 


t 


dt 


If  the  argument  x  +  iy  satisfies  both  of  the  inequalities 

1  <  xz  +  yz  <  1600  —  x  +  0.064  yz>  0 


(39) 


(40) 


then  the  exponential  integral  is  computed  from  the  rational  approximation.  When  the 
converging  factor  is  approximated  by  poles  and  residues,  then  the  exponential  integral 
is  expressed  by  the  equation 

Ei{z)  =  e*  X  (41) 

i=i  z  -  °i 

for  which  the  positions  6i  and  the  residues  of  the  poles  are  listed  in  Table  I. 

If  the  argument  x  +  iy  satisfies  the  inequality 

xz  +  yaSl600  (42) 


then  the  exponential  integral  is  computed  from  the  descending  power  series.  Repeated 
integration  by  parts  leads  to  the  equation 


Ei(z) 


ez  m! 

—  X  -tt+M 

z  __n  z 

771=0 


r 


dt 


(43) 


in  which  the  series  is  asymptotic.  Summation  of  the  series  is  continued  until  there 
is  no  change  in  the  sum  or  until  m  =  39. 

Programming 

SUBROUTINE  CEXPLI  (M0,  AZ,  FE) 

FORTRAN  SUBROUTINE  FOR  COMPLEX  EXPONENTIAL  INTEGRAL 

The  mode  of  operation  is  given  in  M0.  The  real  and  imaginary  parts  of  the  argument 
z  are  given  in  array  AZ.  The  complex  exponential  integral  is  computed  by  series 
expansions  and  rational  approximations.  If  MO  =  0,  the  real  and  imaginary  parts  of 
the  function  Ei(z)  are  stored  in  array  FE.  If  MO  =  1 ,  the  real  and  imaginary  parts  of 
the  function  e~zEi(z)  are  stored  in  array  FE. 
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COMPLEX  FRESNEL  INTEGRAL 


Analysis 

Various  forms  of  complex  integral  with  different  constants  of  integration  can  be 
utilized.  Let  the  complex  Fresnel  integral  E(z )  for  complex  argument  z  be  defined  by 
the  equation 


E(z) 


(44) 


where  the  path  of  integration  lies  within  that  part  of  the  complex  plane  from  which 
the  positive  real  axis  is  excluded.  The  phase  of  z  is  limited  to  the  range  0  to  27r,  and 
the  phase  of  z 1/1,2  is  half  the  phase  of  z. 

The  conventional  Fresnel  integrals  are  defined  in  terms  of  harmonic  functions  by 
the  equations 

C{v)  +  iS(v)  =  f  du  =  f  iu  du  (45) 

Jo  Jo 

or  are  expressed  in  terms  of  Hankel  functions  by  the  equations 

h  f*X  •  f*X 

C(x)  +  iS(x)  =  , —  du  =  “  H±l\u )  du  (46) 

v 27 r  Jo  yZ  ^  Jo  2 


where  x  and  v  are  related  by  the  equation 

1  2 
X  =  g7T^ 


(47) 


The  conventional  integrals  are  recovered  from  the  complex  integral  by  the 
transformations 


t  ->  iu 


z  ->  ix 


and  are  expressed  by  the  equation 


C(x)  +  iS(x)  = 


1  +  i 

-  + 

2 


1  —  i  .  . 
— 7=-  E(ix) 
v  2 


The  transformations 

t  -»  —  uz 

give  the  error  function  of  real  argument 

■•M 


±H(x)  = 


2  r 
V7T  Jo 


e  u  du  =  1  —  i^/2E(-xz) 


where  the  ±  sign  is  the  same  as  the  sign  of  x. 


(48) 


(49) 


(50) 


(51) 
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The  transformations 


(52) 


t  -»  +  uz 


z  -*  +  x 


give  the  error  function  of  imaginary  argument 


T  iH(ix)  = 


V7T 


Jo  ' 


eu  du  =  i  +  ^f2E(+xz) 


where  the  =F  sign  is  opposite  to  the  sign  of  x. 

The  Stieltjes  form  of  the  Fresnel  integral  is  given  by  the  equation 


£(*)  = 


(fz)2e 


:)2e  Cm 

77  Jo  n,k 


du 


uz(z  —  u) 


(53) 


(54) 


both  sides  of  which  approach  zero  as  z  -*  —  and  have  identically  the  same  derivatives 
with  respect  to  z. 

If  the  argument  x  +  iy  satisfies  the  inequality 

xz  +  yz  S  1  (55) 

or  both  of  the  inequalities 

xz  +  yz  <  1444  -  x  +  0.064  yz  £  0  (56) 


then  the  Fresnel  integral  is  computed  from  the  ascending  power  series.  Variation  of 
the  path  of  integration  leads  to  the  equation 


E{z) 


(57) 


Expansion  in  power  series  and  term  by  term  integration  leads  to  the  classical  equation 


i  /2z\a  00  zm 

2?(z)  =  --5=+(— )  D  jz — —r — -  (58) 

V2  \  77  /  m= 0  (2m  +  !)m! 

for  which  the  constant  of  integration  is  given  by  the  gamma  function  T(|)  and  is 
defined  by  the  equation 


e  _ 

— j-  dt  =  VTT 

f2 


(59) 


If  the  argument  x  +  iy  satisfies  both  of  the  inequalities 

1  <  +  yz  <  1444  -  x  +  0.064  yz  >  0  (60) 


then  the  Fresnel  integral  is  computed  from  the  rational  approximation.  When  the 
converging  factor  is  approximated  by  poles  and  residues,  then  the  Fresnel  integral  is 
expressed  by  the  equation 


E(z) 


*  -  <5i 


(61) 


for  which  the  positions  and  the  residues  ei  of  the  poles  are  listed  in  Table  II. 
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If  the  argument  x  +  iy  satisfies  the  inequality 

xz  +  yz>  1444  (62) 


then  the  Fresnel  integral  is  computed  from  the  descending  power  series.  Repeated 
integration  by  parts  leads  to  the  equation 


E(z) 


e*  "-1  (2m)!  (2N)\  (*» 

(2 nzyz  rn=o  22mm!zm  V2n2ZNN\  J_» 


(63) 


in  which  the  series  is  asymptotic.  Summation  of  the  series  is  continued  until  there 
is  no  change  in  the  sum  or  until  m  =  37. 

Programming 

SUBROUTINE  CFRNLI  (MO,  AZ,  FE) 

FORTRAN  SUBROUTINE  FOR  COMPLEX  FRESNEL  INTEGRAL 

The  mode  of  operation  is  given  in  MO.  The  real  and  imaginary  parts  of  the  argument 
z  are  given  in  array  AZ.  The  complex  Fresnel  integral  is  computed  by  series  expansions 
and  rational  approximations.  If  MO  =  0,  the  real  and  imaginary  parts  of  the  function 
E(z)  are  stored  in  array  FE.  If  MO  =  1 ,  the  real  and  imaginary  parts  of  the  function 
e~*E(z)  are  stored  in  array  FE. 
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DOUBLE  EXPONENTIAL  INTEGRAL 


Analysis 

The  double  exponential  integral  Di{z)  is  defined  in  terms  of  its  argument  2  by  the 
equation 


Di(z) 


(64) 


where  the  path  of  integration  lies  in  that  part  of  the  complex  plane  from  which  the 
positive  real  axis  is  excluded. 

If  the  argument  x  +  iy  satisfies  the  inequality 


xz  +  yz  g  1 


(65) 


or  both  of  the  inequalities 

xz  +  yz<  2025  -  x  +  0.064  y2  g  0  (66) 

then  the  double  exponential  integral  is  computed  with  an  ascending  power  series. 
Substitution  of  -t  for  t  and  integration  by  parts  leads  to  the  equation 


r  el 

Di(z)  =  log  (-2)  J  —  dt 
"  (1-0 


log  t  dt  + 


noc 

J  1 


t 


log  t  dt 


log t 
t 


dt  + 


/; 


(1-0 

t 


log  t  dt 


(67) 


for  which  the  constant  of  integration  can  be  evaluated  from  the  limiting  value  of  the 
difference 


tv  le  ‘  log  t  dt 


r  t  v~i 

Jo 


log  t  dt 


as  v  ->  0.  Thus  the  reciprocal  of  v  is  given  by  the  equation 


-=  ft*'-1 

v  Jo 


dt 


and  its  derivative  is  given  by  the  equation 


log  t  dt 


(68) 


(69) 


(70) 
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The  gamma  function  r(y)  is  given  by  the  equation 


'  dt 


(71) 


poo 

r(i/)=  t^e-* 

Jo 

and  its  derivative  is  given  by  the  equation 

poo 

V(v)=  tv~1e~‘  log  t  dt  (72) 

Jo 

From  page  236  of  Reference  23,  the  gamma  function  is  given  also  by  the  equation 


r(u)  = 


n 


v 

1  +  - 
71 


(73) 


and  its  logarithmic  derivative  is  given  by  the  equation 

I»  1  f  1 

T(u)  7  v  l/nt',n(n+  v) 

The  sum  of  the  inverse  squares  of  integers  is  given  by  the  equation 

f  jL  =  zk! 

nz  6 


(74) 


(75) 


Substitution  of  these  formulae  into  the  above  difference  between  integrals  shows  that 
the  limit  as  u  -*  0  is  expressed  by  the  equation 


e-t  pi 

—  log  t  dt  - 
o  £  Jo 


log  t  rr2  72 

- dt  =  — -  H - 

t  12  2 


Expansion  in  series  and  term  by  term  integration  leads  to  the  equation 

«(,)^  +  |(7  + lo,  (-.)).  ^-^ 

in  which  the  series  is  absolutely  convergent. 

If  the  argument  x  +  iy  satisfies  both  of  the  inequalities 


(76) 


(77) 


1  <xz  +  yz<  2025 


x  +  0.064  yz  >  0 


(78) 


then  the  double  exponential  integral  is  computed  from  the  rational  approximation. 
When  the  converging  factor  is  approximated  by  poles  and  residues,  then  the  double 
exponential  integral  is  expressed  by  the  equation 


Di(z)  =  —  E 


z  i=1  z  -  <5i 

for  which  the  positions  <5t  and  the  residues  e4  of  the  poles  are  listed  in  Table  III. 


(79) 
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If  the  argument  x  +  iy  satisfies  the  inequality 

xz  +  yz>  2025 


(80) 


then  the  double  exponential  integral  is  computed  from  the  descending  power  series. 
Expansion  of  the  integrand  as  a  descending  series  and  progressive  integration  by  parts 
gives  the  approximation 


Di(z) 


N-l 


£ 


(m+  1)!  “ 

zm  *.„  *  +  1 


(81) 


in  which  the  series  is  asymptotic.  Summation  of  the  series  is  continued  until  there 
is  no  change  in  the  sum  or  until  m  =  44. 

Programming 

SUBROUTINE  CDEXPI  (MO,  AZ,  FD) 

********************************************************************************************* 
FORTRAN  SUBROUTINE  FOR  DOUBLE  EXPONENTIAL  INTEGRAL 

**5|<**:t:**5fc*************************************5)i:************:+:*********:t:**:t::|c*****jtc*:4;**:5t::|<::+:****** 


The  mode  of  operation  is  given  in  MO.  The  real  and  imaginary  parts  of  the  argument 
z  are  given  in  array  AZ.  The  double  exponential  integral  is  computed  by  series  expansions 
and  rational  approximations.  If  MO  =  0,  the  real  and  imaginary  parts  of  the  function 
Di(z)  are  stored  in  array  FD.  If  MO  =  1,  the  real  and  imaginary  parts  of  the  function 
e~zDi(z)  are  stored  in  array  FD. 
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DOUBLE  FRESNEL  INTEGRAL 


Analysis 

The  double  Fresnel  integral  D(z)  is  defined  in  terms  of  its  argument  z  by  the 
equation 


D(z) 


(82) 


where  the  path  of  integration  lies  within  that  part  of  the  complex  plane  from  which 
the  positive  real  axis  is  excluded.  The  phase  of  z  is  limited  to  the  range  0  to  2it,  and 
the  phase  of  z1^2  is  half  the  phase  of  z. 

If  the  argument  x  +  iy  satisfies  the  inequality 


x2  +  y2<  1 


(83) 


or  both  of  the  inequalities 


x2  +  y2  <  1849  -  x  +  0.064  y2  ^  0  (84) 

then  the  double  Fresnel  integral  is  computed  from  the  ascending  power  series. 
Substitution  of  —  t  for  t  and  integration  by  parts  leads  to  the  equation 

r  e‘ 

D(z)  =  log  (-z)  -j  dt 
J-=o  ^2 


+  i 


-  i 


e 

— j-  log  t  dt 
1 2 

:  e~l 

— j-  log  t  dt 
t ® 


(85) 


for  which  the  constant  of  integration  is  given  by  the  equation 

n 


"'(f)  =  f  ^~T  log  i  dt  ~  -  (7  +  2  log  2)Vrr 
Jo  tz 


(86) 


Expansion  in  series  and  term  by  term  integration  leads  to  the  equation 


D(z)  =  —  i( 7  +  2  log  2  +  log  (— z))V7T  +  Y, 

m=0 


(m  +  |)2m! 


(8.7) 


in  which  the  series  is  absolutely  convergent. 
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If  the  argument  x  +  iy  satisfies  both  of  the  inequalities 


1  <  +  yz  <  1849  -  x  +  0.064  y2  >  0  (88) 

then  the  double  Fresnel  integral  is  computed  from  the  rational  approximation.  When 
the  converging  factor  is  approximated  by  poles  and  residues,  then  the  double  Fresnel 
integral  is  expressed  by  the  equation 

ez  19  e. 

(89) 

Z*  i=l  *  -  <*i 


for  which  the  positions  <5*  and  the  residues  et  of  the  poles  are  listed  in  Table  IV. 

If  the  argument  x  +  iy  satisfies  the  inequality 

xz  +  yz>  1849  (90) 


then  the  double  Fresnel  integral  is  computed  from  the  descending  power  series. 
Expansion  of  the  integrand  as  a  descending  series  and  progressive  integration  by  parts 
gives  the  approximation 


D(z) 


ez_  "-1  (2m  +  1)!  ™  1 

zf  m=o  22mm!zm  k%  2k  +  1 


(91) 


in  which  the  series  is  asymptotic.  Summation  of  the  series  is  continued  until  there 
is  no  change  in  the  sum  or  until  m  =  42. 

Programming 

SUBROUTINE  CDFRNI  (MO,  AZ,  FD) 

FORTRAN  SUBROUTINE  FOR  DOUBLE  FRESNEL  INTEGRAL 

The  mode  of  operation  is  given  in  MO.  The  real  and  imaginary  parts  of  the  argument 
z  are  given  in  array  AZ.  The  double  Fresnel  integral  is  computed  by  series  expansions 
and  rational  approximations.  If  MO  =  0,  the  real  and  imaginary  parts  of  the  function 
D(z)  are  stored  in  array  FD.  If  MO  =  1,  the  real  and  imaginary  parts  of  the  function 
e~zD(z)  are  stored  in  array  FD. 
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EXPONENTIAL  EXPONENTIAL  INTEGRAL 


Analysis 

The  exponential  exponential  integral  k(z)  is  defined  in  terms  of  its  argument  z  by 
the  equation 


e~“  pu  ei 

-  —  dt  du  (92) 

U  J-o o  t 

where  the  path  of  integration  lies  in  that  part  of  the  complex  plane  from  which  the 
positive  real  axis  is  excluded.  Insofar  as  the  exponential  integral  approaches  its 
asymptotic  approximation  for  large  argument,  the  integrand  and  its  integral  approach 
zero  along  any  circle  of  infinite  radius  which  is  centered  at  the  origin.  When  the  path 
of  integration  is  deformed  to  follow  partly  an  arc  of  infinite  radius,  the  integration 
from  —  °°  to  2  is  converted  into  the  negative  of  the  integration  from  z  to  +°°.  The  path 
of  integration  must  not  cross  the  positive  real  axis. 

If  the  argument  x  +  iy  satisfies  the  inequality 


x2  +  yz  <  1 


(93) 


then  the  exponential  exponential  integral  is  given  by  the  equation 


k(z) 

=  irr2 

+  1(7  +  log  ( ~z))z 

OO 

+  E 

^7  +  log  (- 

m=  1 

mf  mm\ 

OO 

-  E 

(~z)m  £ 

1 

m=  1 

mm!  k=i 

k 

in  which  the  series  is  convergent. 

If  the  argument  x  +  iy  satisfies 

both 

of  the  inequalities 

1  <x2  +  y2<  1296 

-  x  +  0.070 

then  the  exponential  exponential  integral  is  given  by  the  equation 


(94) 


(95) 


k{z)  =  (7  +  log  (—z))Ei(—z) 


00  <»  -1 

-Di(-z)-e~*  £  —  E  ~z 

m= 0  m'  k=m+ 1 


(96) 


for  which  the  integrals  Ei(—z)  and  Di(—z)  are  computed  from  their  rational 
approximations.  Of  two  alternative  series  expansions  the  one  given  above  has  the 
advantage  that  the  terms  are  all  of  the  same  sign  on  the  positive  real  axis. 
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If  the  argument  x  +  iy  satisfies  both  of  the  inequalities 

1  <xz  +  yz<  1296  -  x  +  0.070  y2  >  0 

then  the  exponential  exponential  integral  is  given  by  the  equation 

18 

k(z)  =  -  £ 


z  —  <5,- 


(9?) 


(98) 


for  which  the  positions  <5*  and  the  residues  ei  of  the  poles  are  listed  in  Table  V. 
If  the  argument  x  +  iy  satisfies  the  inequality 


xz  +  yz>  1296 


(99) 


then  the  exponential  exponential  integral  is  given  by  the  approximation 


k(z) 


<  AT— 1 
-  -  £ 


Z 


m- 0 


m! 

(m  +  1  )zm 


(100) 


in  which  the  series  is  asymptotic.  Summation  of  the  series  is  continued  until  there 
is  no  change  in  the  sum  or  until  m  =  35. 

Pro  gramming 

SUBROUTINE  CEXEXI  (AZ,  FK) 

***********************+:**+***+**+%******************+*************************************** 
FORTRAN  SUBROUTINE  FOR  EXPONENTIAL  EXPONENTIAL  INTEGRAL 


The  real  and  imaginary  parts  of  the  argument  z  are  given  in  array  AZ.  The  exponential 
exponential  integral  is  computed  by  series  expansions  and  rational  approximations. 
The  real  and  imaginary  parts  of  the  function  k{z)  are  stored  in  array  FK. 
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EXPONENTIAL  FRESNEL  INTEGRAL 


Analysis 

The  exponential  Fresnel  integral  m(z )  is  defined  in  terms  of  its  argument  z  by  the 
equation 


rx  e~u  Cu  el 

m(z)  =  — y  — dt  du  ( 1 0 1 ) 

tl— OO  */ —  oo  ^ 

where  the  path  of  integration  lies  in  that  part  of  the  complex  plane  from  which  the 
positive  real  axis  is  excluded.  The  phase  of  z  is  limited  to  the  range  0  to  2ix,  and  the 
phase  of  z1/2  is  half  the  phase  of  z.  Insofar  as  the  exponential  integral  approaches 
its  asymptotic  approximation  for  large  argument,  the  integrand  and  its  integral 
approach  zero  along  any  circle  of  infinite  radius  which  is  centered  at  the  origin.  When 
the  path  of  integration  is  deformed  to  follow  partly  an  arc  of  infinite  radius,  the 
integration  from  —  °°  to  z  is  converted  into  the  negative  of  the  integration  from  z  to 
+  °°.  The  path  of  integration  must  not  cross  the  positive  real  axis. 

If  the  argument  x  +  iy  satisfies  the  inequality 

xz+yz£l  (102) 

then  the  exponential  Fresnel  integral  is  given  by  the  equation 

3 

m(z)  =  in2 


+  za  £  (j  +  log  (-z)  -  -71 

m=0  \  m+Z 


1  \  (~zY 


(m  +  |)m! 


_  I  y  (~*)m  y  L 
rti  (m  +  |)m!  k=l  k 

in  which  the  series  is  convergent. 

If  the  argument  x  +  iy  satisfies  both  of  the  inequalities 

1  <  xz  +  yz  <  1296  -  x  +  0.070  y2  g  0 

then  the  exponential  Fresnel  integral  is  given  by  the  equation 


(103) 


(104) 


m(z)  =  =F  (7  +  2  log  2  +  log  (— z)) y/2niE(—z) 


±  iD(-z)  -  e  *  £ 

m= 0 


m+A 

*5 


r(m  + 1) 


£ 

ifc=m+l 


r(fc  +  |) 

kk\ 


(105) 


for  which  the  integrals  E{—z)  and  D{-z)  are  computed  from  their  rational 
approximations.  The  upper  sign  is  valid  when  y  >  0  and  the  lower  sign  is  valid  when 
y  ^  0.  Of  two  alternative  series  expansions  the  one  given  above  has  the  advantage  that 
the  terms  are  all  of  the  same  sign  on  the  positive  real  axis. 
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If  the  argument  x  +  iy  satisfies  both  of  the  inequalities 


1  <xz  +  yz<  1296  -  x  +  0.070  yz  >  0 

then  the  exponential  Fresnel  integral  is  given  by  the  equation 

i  is  *. 

•771(2)  =  -  2z2  E 


(106) 


5, 


(107) 


for  which  the  positions  8i  and  the  residues  ei  of  the  poles  are  listed  in  Table  VI. 
If  the  argument  x  +  iy  satisfies  the  inequality 


xz  +  yz>  1296 

then  the  exponential  Fresnel  integral  is  given  by  the  approximation 


771(2) - 1  E 


1  N~l  m\ 


I  m=o  (m  +  kW 


(108) 


(109) 


in  which  the  series  is  asymptotic.  Summation  of  the  series  is  continued  until  there 
is  no  change  in  the  sum  or  until  m  =  35. 

Programming 

SUBROUTINE  CEXFRI  (AZ,  FM) 

FORTRAN  SUBROUTINE  FOR  EXPONENTIAL  FRESNEL  INTEGRAL 

The  real  and  imaginary  parts  of  the  argument  2  are  given  in  array  AZ.  The  exponential 
Fresnel  integral  is  computed  by  series  expansions  and  rational  approximations.  The 
real  and  imaginary  parts  of  the  function  m(z)  are  stored  in  array  FM. 
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LOGARITHMIC  EXPONENTIAL  INTEGRAL 


Analysis 

The  logarithmic  exponential  integral  K(z,  q )  is  defined  by  the  equation 


K(z,  q) 


fj 


log  1  H —  )  dt 

q/ 


(no) 


in  which  z  is  the  argument  and  q  is  a  parameter.  The  integrand  is  analytic  at  t  —  0 
but  has  a  singularity  at  t  —  —q.  The  complex  plane  is  cut  along  a  line  which  extends 
from  the  singularity  to  infinite  distance  in  a  radial  direction. 

The  logarithm  for  complex  argument  may  be  defined  by  the  equation 


log  1  + 


t  \  r‘  du 
qj  J0  q  +  u 


where  the  path  of  integration  extends  radially  outward  from  the  origin. 
If  the  argument  x  +  iy  satisfies  the  inequality 


xz  +  yz  <  1 


or  both  of  the  inequalities 
xz  +  yz<  1600 


—  x  +  0.064  yz  £  0 


(111) 


(112) 


(113) 


then  the  logarithmic  exponential  integral  is  computed  from  an  ascending  power  series. 
The  path  of  integration  may  be  varied  to  pass  through  the  origin,  and  the  exponential 
function  in  the  integrand  may  be  expanded  in  a  power  series. 

The  constant  of  integration  is  given  by  the  equation 


e  /  t 

—  log  1  +  - 
t  \  q 


dt  =  — 


dt  du 


(114) 


both  sides  of  which  approach  zero  as  q  -»  and  have  identically  the  same  derivatives 
with  respect  to  q. 

The  logarithmic  exponential  integral  is  given  by  the  equation 


K(z,q)  =  -k(q)+L  1  +  -  +  E  T„ 

7/  m=l 


for  which  the  terms  are  defined  by  the  equation 


i  (  „  /  z\  r  um 

Tm  = - r)zmlog  1+-  -  — —du 

m  m\  {  \  q  J  J0  q  +  u 


(115) 


(116) 
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When  |g|  s>  V1.3|z|  the  series  is  started  with  the  equation 


Ti  =  (q  +  z)  log  1  +  -  -  0 
9/ 


(117) 


and  is  continued  in  the  direction  of  ascending  order  with  the  recurrence  equation 

^m—1  / 


T  = 

A  m 


){q  +  z)  log 


H) 


m , 


(m  -  1) 

2  9^m- 1 


(118) 


The  recurrence  is  cycled  to  the  order  N  to  which  the  ascending  series  of  the  exponential 
integral  must  be  carried  for  accurate  evaluation. 

When  |qr]  >  VI. 3 \z\  the  terms  are  given  by  the  equation 


m  °° 

—x  E 

m-  k=i 


i 

k(m  +  k) 


(119); 


The  series  is  started  at  the  order  N  and  is  continued  in  the  direction  of  descending 
order  with  the  recurrence  equation 


T  = 

■*  m 


mm! 


1  2 
m  +  1  q 


(m  +  l)2 
mg 


^m+l 


(120) 


The  terms  in  the  expansion  of  the  exponential  integral  are  computed  in  advance  of 
recurrence  and  are  stored  in  an  array. 

An  integration  by  parts  leads  to  the  equation 


dt  =  log 


du 
q  +  u 


(121) 


If  the  modulus  of  z  is  large  enough,  the  exponential  integral  can  be  replaced  by  its 
rational  or  asymptotic  approximations.  The  path  of  integration  must  lie  outside  of  a 
circle  of  radius  \z\.  If  this  path  would  cross  a  cut,  then  the  actual  integrand  must  be 
integrated  in  along  the  cut,  around  the  singularity,  and  out  along  the  cut.  To  the 
integration  along  the  outer  path  must  be  added  the  correction 

r  e*  rg  e‘  ,  , 

±  Zm  \  —  dt  T  Zm  I  -  di  (122) 


where  the  upper  sign  is  used  when  the  singularity  is  above  the  real  axis  and  the  lower 
sign  is  used  when  the  singularity  is  on  or  below  the  real  axis. 

In  the  integration  along  the  outer  path  the  imaginary  part  of  the  logarithm  does 
not  reverse  sign  where  the  path  crosses  the  cut.  If  the  logarithm  is  given  its  usual 
value  at  the  end  of  integration,  then  an  additional  correction  in  the  amount  of 


T  27Ti 


(123) 


must  be  added  to  the  integration  along  the  outer  path.  Only  the  integration  to  —  q 
remains  when  both  corrections  are  added. 
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A  crossing  of  the  cut  occurs  only  if  \z\  >  |g|  and  the  inequalities 


3m  q<  0  3m  z>  0  3m  log  ^1  +  —  j  >  0  (124) 

or  the  inequalities 

3m  q>  0  3mz1k  0  ,7m  log  ^lH — j  <  0  (125) 

also  are  satisfied.  If  q  happens  to  lie  on  the  positive  real  axis,  then  a  crossing  occurs 
if  the  inequality 


3m  log 


g  0 


(126) 


is  satisfied. 

If  the  argument  x  +  iy  satisfies  both  of  the  inequalities 

1  <  xz  +  yz  <  1600  -  x  +  0.064  yz  >  0 


(127) 


then  the  logarithmic  exponential  integral  is  computed  from  the  rational  approximation 
for  the  exponential  integral.  Along  the  outer  path  the  contribution  to  the  integration 
is 


K{z,  q)  ~  £  et 

t-1 


e  /  z  \  e  x 

Z  -  <5i  l0g  \ 1  +  q  )  q  +  (5 


"*2 

i  v  - 


z~ei  e‘ 


dt  H -  , 

t  g  +  <5i  J-o. 


dt 


(128) 


for  which  the  positions  6i  and  the  residues  e{  of  the  poles  are  listed  in  Table  I. 

If  the  outer  path  crosses  the  cut,  and  also  ends  closer  to  the  real  axis  than  the 
singularity,  or  if  the  outer  path  does  not  cross  the  cut,  yet  ends  farther  from  the 
real  axis  than  the  singularity,  then  a  correction  must  be  added  to  that  exponential 
integral  whose  argument  is  q  +  z.  The  correction  is 

±2t tx  (129) 


where  the  sign  is  determined  by  the  direction  in  which  the  singularity  must  be  encircled 
in  order  to  correct  the  evaluation  of  the  exponential  integral.  If  \z\  >  |q|,  then  the  lower 
sign  is  used  when  the  singularity  is  above  the  real  axis  and  3m(q  +  z)  g  0,  but  the 
upper  sign  is  used  when  the  singularity  is  below  the  real  axis  and  3m{q  +  z)  >  0.  If 
\z\  ^  |g|,  then  the  signs  are  reversed. 

If  the  arguments  satisfy  both  of  the  inequalities 

|g  +  <5t|gf  |g  +  <5t|  ^  f|g  +  z|  (130) 


then  the  last  two  integrals  in  the  rational  approximation  are  nearly  equal.  Their 
difference  can  be  expressed  by  a  series  expansion  in  accordance  with  the  equation 


1 


q  +  L 


-<*-«*> 


Z-6.  t  P9  +  *  r,t 

-  dt  -  - 

J  —  CO  t 


t 


dt 


=  Y,Tn 


(131) 
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where  the  term  of  lowest  order  is  given  by  the  equation 


To  =  ~ 


q  +  z 


+  e-(q+z) 


r+z  - 

J—  t 


dt 


and  the  terms  of  higher  order  are  given  by  the  equations 

(O'  +  O 


T  1  (q  +  ( q  +  gil  t 

71  (n  +  l)(q  +  z)  \  q  +  z  /  (n  +  1)  n_1 


The  recurrence  is  cycled  until  there  is  no  change  in  the  summation. 
If  the  argument  x  +  iy  satisfies  the  inequality 

xz  +  yz>  1600 


(132) 


(133) 


(134) 


then  the  logarithmic  exponential  integral  is  computed  from  the  asymptotic 
approximation  for  the  exponential  integral.  Along  the  outer  path  the  contribution  to 
the  integration  is 


n-  1 

K(z,  q)  ~  ez  E  Tn 

71=  1 


(135) 


where  the  terms  are  defined  by  the  equation 
(n  -  1)! 


T  = 

4  T). 


l0g(l  +  i)_(„_1)!e-.J___£l_,„  (136) 


When  \z\  £  V6|g|  the  series  is  started  with  the  equation 


T,  =  -  log  (l  +  -)  -  -  f  —dt+~ 
z  \  q)  q  t  q 


,-(9+3) 


t 


dt 


(137) 


and  is  continued  in  the  direction  of  descending  order  with  the  recurrence  equation 


T  = 

4  n. 


(n-  1)! 


i  +  -)iog(i  +  -)-h— tl!e- 

q!  \  q)  q 


r  el 
tn 


(n  -  1) 


Tn-!  (138) 


When  \z\  >  V6|q|  the  series  is  started  at  the  order  N  and  is  continued  in  the  direction 
of  ascending  order  with  the  recurrence  equation 


T  = 

4  Tt. 


(n-  1)! 


1)!  .  .  /  z\  .  .  e*  q 

T~  (q  +  z)  log  ( 1  +  —  )  —  (n  -  1)!  e  — n  dt  -  -  Tn+1 

\  q  /  j -00  ®  n 


(139) 


In  either  case  the  exponential  integral  is  generated  with  the  recurrence  equation 

^  (n  —  1)! 


(n  —  l)!e' 


f  e‘ 
J-oc  t 


+  n\e 


r  jL 


dt 


(140) 


The  terms  for  the  asymptotic  expansion  are  computed  in  advance  of  recurrence  and 
are  stored  in  an  array.  The  maximum  order  N  is  determined  by  whichever  of  the 
criteria 


N\ 


«  1 


<<  1 


(141) 
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is  appropriate  for  accurate  computation. 

When  q  +  z  =  0,  the  logarithmic  exponential  integral  is  taken  to  be  that  value  which 
is  approached  when  z/q  is  real  and  approaches  —  1. 

Programming 

SUBROUTINE  CLGEXI  (MO,  AZ,  CQ,  FK) 

********************************************************************************************* 
FORTRAN  SUBROUTINE  FOR  LOGARITHMIC  EXPONENTIAL  INTEGRAL 
********************************************************************************************* 

The  mode  of  operation  is  given  in  MO.  The  real  and  imaginary  parts  of  the  argument 
2  are  given  in  array  AZ.  The  real  and  imaginary  parts  of  the  parameter  q  are  given 
in  array  CQ.  The  logarithmic  exponential  integral  is  computed  by  series  expansions 
and  rational  approximations.  Calls  are  made  to  Subroutine  CLGMCI,  Subroutine  CEXPLI, 
and  Subroutine  CEXEXI.  If  MO  =  0,  the  real  and  imaginary  parts  of  the  function  K(z,  q) 
are  stored  in  array  FK.  If  MO  =  1,  the  real  and  imaginary  parts  of  the  function  e~zK(z ,  q) 
are  stored  in  array  FK. 
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ARCTAN GENTIAL  FRESNEL  INTEGRAL 


Analysis 

The  arctangential  Fresnel  integral  M(z,  q)  is  defined  by  the  equation 


(142) 


in  which  z  is  the  argument  and  q  is  a  parameter.  The  integrand  is  analytic  at  t  =  0 
but  has  a  singularity  at  t  =  —q.  The  complex  plane  is  cut  along  a  line  which  extends 
from  the  singularity  to  infinite  distance  in  a  radial  direction. 

The  arctangent  for  complex  argument  may  be  defined  by  the  equation 


tan  1 


du 

~~T 

uz(q  +  u ) 


(143) 


where  the  path  of  integration  extends  radially  outward  from  the  origin.  Resolution  of 
the  integrand  into  partial  fractions  and  integration  leads  to  the  equation 

tan_1(g)  =  ^log  I1  -i(“)  J  j1  +i(“)  J  (144) 

which  is  useful  for  numerical  evaluation. 

If  the  argument  x  +  iy  satisfies  the  inequality 

xz  +  yz£  1  (145) 

or  both  of  the  inequalities 

xz  +  yz<  1444  -a:  +  0.064y8g0  (146) 


then  the  arctangential  Fresnel  integral  is  computed  from  an  ascending  power  series. 
The  path  of  integration  may  be  varied  to  pass  through  the  origin,  and  the  exponential 
function  in  the  integrand  may  be  expanded  in  a  power  series. 

The  constant  of  integration  is  given  by  the  equation 


e 

—  tan 
tz 


i 


dt=- 


—  dt  du 


d4?) 


both  sides  of  which  approach  zero  as  q  -*  —  °°  and  have  identically  the  same  derivatives 
with  respect  to  q. 


27 


The  arctangential  Fresnel  integral  is  given  by  the  equation 

oo 

M(z,  q)  =  -  | m(q)  +  E  L 


for  which  the  terms  are  defined  by  the  equation 


"  (m  +  iJmlf  W  ”  Jo  9 

When  |g|  ^  V1.3|z|  the  series  is  started  with  the  equation 


(i g2  -  izz)  log  \  1  -  i\ 


(g2  +  izz)  log  ( 1  +  i 


and  is  continued  in  the  direction  of  ascending  order  with  the  recurrence  equation 


(m  +  |)m! 


(g  +  z) tan 


JLf£ 

V  g  /  2m\g 


I)  t 

(m  +  |)  m  m_1 


The  recurrence  is  cycled  to  the  order  N  to  which  the  ascending  series  of  the  Fresnel 
integral  must  be  carried  for  accurate  evaluation. 

When  |g|  >  VI.3|z|  the  terms  are  given  by  the  equation 

1  zm  °°  1  /  z\k 

Tm=~ ^  R  (fc-D(m^)h)  (152) 

The  series  is  started  at  the  order  N  and  is  continued  in  the  direction  of  descending 
order  with  the  recurrence  equation 


(m  +  f)m! 


1  +  —  1  tan  1  —  - 


m  +  1  \  g 


(m  +  |)  (m  +  1) 
(m  +  |)  g 


Tm+1  (153) 


The  terms  in  the  expansion  of  the  Fresnel  integral  are  computed  in  advance  of 
recurrence  and  are  stored  in  an  array. 

An  integration  by  parts  leads  to  the  equation 


f  ^  tan-^y  df  =  tan-(if  f  ij  *  -  fr*  f 
J-ooA  \<?/  W/  J-co  A  J- oo0,I 


^  dt  (154) 


tz  w/  w/  w2(g  +  w)  f2 

If  the  modulus  of  z  is  large  enough,  the  Fresnel  integral  can  be  replaced  by  its  rational 
or  asymptotic  approximations.  The  path  of  integration  must  lie  outside  of  a  circle  of 
radius  |z|.  If  this  path  would  cross  a  cut,  then  the  actual  integrand  must  be  integrated 
in  along  the  cut,  around  the  singularity,  and  out  along  the  cut.  To  the  integration 
along  the  outer  path  must  be  added  the  correction 


-*  e*  r9  e‘ 

7T  —  dt  ±  TT  I  -j 

U  —co  J.  2  J —oo  4-  2 


where  the  lower  sign  is  used  when  the  singularity  is  on  the  negative  real  axis,  but  the 
upper  sign  is  used  when  the  singularity  is  anywhere  else. 
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In  the  integration  along  the  outer  path  the  real  part  of  the  arctangent  does  not 
reverse  sign  where  the  path  crosses  the  cut.  If  the  arctangent  is  given  its  usual  value 
at  the  end  of  integration,  then  an  additional  correction  in  the  amount  of 


(156) 


must  be  added  to  the  integration  along  the  outer  path.  Only  the  integration  to  —q 
remains  when  both  corrections  are  added. 

A  crossing  of  the  cut  occurs  only  if  \z\  >  |gj  and  the  inequalities 

1 

Smq<  0  9mz>  0  tan-1^— j  <0  (157) 

or  the  inequalities 

i 

9m  q  >  0  3mz^  0  %e  tan-1^— ^  <0  (158) 

also  are  satisfied.  If  q  happens  to  lie  on  the  positive  real  axis,  then  a  crossing  occurs 
if  the  inequality 


%e  tan  — j 


&  0 


(159) 


is  satisfied. 

If  the  argument  x  +  iy  satisfies  both  of  the  inequalities 


1  <xz  +  yz<  1444  -x  +  0.064  yz  >  0  (160) 

then  the  arctangential  Fresnel  integral  is  computed  from  the  rational  approximation 
for  the  Fresnel  integral.  Along  the  outer  path  the  contribution  to  the  integration  is 


M(z,  q )  ■ 


18 

E 

i=  1 


<5,- 


i 

tan-1!  —  )  — 


i  \  «-■ 
z^e 


:  p  il  +  b!£!  f 

i  J- co  t  q  +  6t  J. 


q  +  6. 


—  dt 
t 


(161) 


for  which  the  positions  and  the  residues  ei  of  the  poles  are  listed  in  Table  II. 

If  the  outer  path  crosses  the  cut,  and  also  ends  closer  to  the  real  axis  than  the 
singularity,  or  if  the  outer  path  does  not  cross  the  cut,  yet  ends  farther  from  the 
real  axis  than  the  singularity,  then  a  correction  must  be  added  to  that  exponential 
integral  whose  argument  is  q  +  z.  The  correction  is 

±2ni  (162) 


where  the  sign  is  determined  by  the  direction  in  which  the  singularity  must  be  encircled 
in  order  to  correct  the  evaluation  of  the  exponential  integral.  If  |z|  >  |g|,  then  the  lower 
sign  is  used  when  the  singularity  is  above  the  real  axis  and  9m(q  +  z)  <  0,  but  the 
upper  sign  is  used  when  the  singularity  is  below  the  real  axis  and  9m(q  +  z)  >  0.  If 
\z\  £  |g|,  then  the  signs  are  reversed. 
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If  the  arguments  satisfy  both  of  the  inequalities 


\g  +  <5il 


\q  +  <5il  ^  fig  +  z| 


then  the  last  two  integrals  in  the  rational  approximation  are  nearly  equal.  Their 
difference  can  be  expressed  by  a  series  expansion  in  accordance  with  the  equation 


g  + 


r*“4t  e‘ 


dt  —  e" 


dt  =  E  Tn 


where  the  term  of  lowest  order  is  given  by  the  equation 


and  the  terms  of  higher  order  are  given  by  the  equations 


1  /  g  +  (h\n  (g  +  (b) 

(n  +  l)(g  +  z)  \  g  +  z  /  +  (n  +  1)  Tn~l 


The  recurrence  is  cycled  until  there  is  no  change  in  the  summation. 

If  the  argument  x  +  iy  satisfies  the  inequality 

x2  +  yz>  1444  (167) 

then  the  arctangential  Fresnel  integral  is  computed  from  the  asymptotic  approximation 
for  the  Fresnel  integral.  Along  the  outer  path  the  contribution  to  the  integration  is 

N- 1 

M(z,  q)~ez  E  Tn  (168) 


where  the  terms  are  defined  by  the  equation 


Tn  ~  tan-1(— 

r(f)zn~* 


rfo-f)  i 


(g  +  u) 


When  \z\  ^  V6|g|  the  series  is  started  with  the  equation 


T,  -  4  )*  -  4 r  J 

~Z  W/  „Z  J -OO  t 


dt+  \  e~(g+z)  —  dt 

qz  J-o o  t 


and  is  continued  in  the  direction  of  descending  order  with  the  recurrence  equations 


r (n  -  i) 


1\  1  nz 

_ ZJ_  2 

i)  i e 

2  /  nZ  J 


(n-I) 


Tn- 1  (171) 


'2/  qZ 
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When  \z\  >  V6  |g|  the  series  is  started  at  the  order  N  and  is  continued  in  the  direction 
of  ascending  order  with  the  recurrence  equation 


r(n-j) 

r(|)2n+l 


(q  +  z)  tan' 


if.. 


r(n  -|).l 


/: 


dt  — 


n  — 


(172) 


In  either  case  the  exponential  integral  is  generated  with  the  recurrence  equation 


r(n-|) 

r(!) 


r(n  - 1) 

r(|)2n 


+ 


(™- f) 


r(n  + 1) 
r(|)  e 


(173) 


The  terms  for  the  asymptotic  expansion  are  computed  in  advance  of  recurrence  and 
are  stored  in  an  array.  The  maximum  order  N  is  determined  by  whichever  of  the 
criteria 


W  + 1) 

r(|)|zr 


<<  1 


<<  i 


(174) 


is  appropriate  for  accurate  computation. 

When  q  +  z  =  0,  the  arctangential  Fresnel  integral  is  taken  to  be  that  value  which 
is  approached  when  z/q  is  real  and  approaches  —1. 

Programming 

SUBROUTINE  CATFRI  (MO,  AZ,  CQ,  FM) 

********************************************************************************************* 
FORTRAN  SUBROUTINE  FOR  ARCTANGENTIAL  FRESNEL  INTEGRAL 
********************************************************************************************* 


The  mode  of  operation  is  given  in  MO.  The  real  and  imaginary  parts  of  the  argument 
z  are  given  in  array  AZ.  The  real  and  imaginary  parts  of  the  parameter  q  are  given 
in  array  CQ.  The  arctangential  Fresnel  integral  is  computed  by  series  expansions  and 
rational  approximations.  Calls  are  made  to  Subroutine  CEXPLI,  Subroutine  CFRNLI,  and 
Subroutine  CEXFRI.  If  MO  =  0,  the  real  and  imaginary  parts  of  the  function  M(z,  q)  are 
stored  in  array  FM.  If  MO  =  1,  the  real  and  imaginary  parts  of  the  function  e~xM(z,q) 
are  stored  in  array  FM. 
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ORDINARY  BESSEL  FUNCTION 


Analysis 

The  Bessel  function  Jv(z)  is  given  by  the  equation 


Jv(z)  =  E 


1  ^\v+2m 


m! r(y  +  m  +  1) 

while  the  Weber  function  Yv(z)  is  given  by  the  equation 

„  ,  N  Jv(z)  cos  vn  -  J-V(z)  . 

K(z)  = - : - - 

sm  UTT 

The  largest  term  in  the  ascending  series  for  Jv(z)  occurs  where 

hM2 


m\v  +  m\ 

The  magnitude  of  the  largest  term  is  given  by  the  Stirling  approximation 

V 

1  i  m 


\h\u+Zm 


m 


!r(i/  +  m+l)  2rrVm(i/  +  m)  \^  + 


m 


,v+2m 


(175) 


(176) 


(177) 


(178) 


which  remains  bounded  as  v  ->  »  if  the  ratio  between  m  and  v  does  not  exceed  a 
maximum  limit.  The  maximum  ratio  a  satisfies  the  transcendental  equation 


a 


1  +  a 


,l+2a  _ 


=  1 


(179) 


whose  solution  is 


a  =  9.983932012886692x10" 


(180) 


The  order  and  the  argument  are  bounded  by  the  limiting  relations 

m  =  av  \z\  =  yv  (181) 

for  which  the  ratio  \x  is  given  by  the  equation 

i 

H  =  2  (a  +  a2)2  =  0.6627434193491816  (182) 

If  v  is  replaced  by  - v ,  then  a  is  replaced  by  -  1  -  a,  while  /z  remains  the  same.  Thus 
the  maximum  term  remains  bounded  if  the  order  and  the  argument  satisfy  the 
inequality 


\z\<y\v\ 


(183) 
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A  popular  method  for  the  computation  of  Bessel  functions  in  the  critical  region 
which  cannot  be  reached  by  the  ascending  series  utilizes  the  recurrence  equations 


2v  .  . 

—  Jv(z) 

(184) 

z 

Zv 

—  Yv{z) 
z 

(185) 

which  are  cycled  in  the  direction  of  descending  order.  During  the  \x th  cycle  a  rounding 
error  is  introduced  and  this  rounding  error  then  persists  in  subsequent  cycles.  In 
view  of  the  recurrence  relation 

Jv(z)Yv+l(z)  -  Jv+i(z)Yv(z)  =  -  —  (186) 

7T2: 

the  persisting  error  in  the  uth  cycle  is  given  by  the  expression 

txz  (  1 

-~Y^{z)Jv{z)-J^+l{z)Yv{z){)Etl  (187) 

for  a  descending  recurrence.  Interchange  of  n  and  v  in  the  coefficient  of  the  expression 
gives  the  persisting  error  for  an  ascending  recurrence.  The  recurrence  can  be  cycled 
in  either  direction  without  decay  or  growth  of  error  wherever  |/„(z)[  and  1 3^, ( z ) |  both 
are  bounded.  The  practical  range  for  recurrence  has  been  determined  by  actual 
computation  to  be  where  \z\  >  0.8|i/|.  Elsewhere  the  recurrence  tends  to  accentuate 
whichever  of  the  functions  is  increasing  relative  to  the  other.  A  large  number  of  cycles 
is  required  when  the  order  is  small  but  the  argument  is  large. 

If  the  order  and  the  argument  satisfy  the  inequality 

|z|  17.5  +  fM2  (188) 

then  the  Bessel  function  is  given  by  the  equation 


Jv(z)  =  ( -  )  }Pv(z)  cos  (z  -  fl/7T  -  \tt)  -  Qv(z)  sin  (z  -  gl/7T  -  j7t)( 


where  the  function  Pv(z)  is  the  sum  of  the  even-ordered  terms  and  the  function  f?„(z) 
is  the  sum  of  the  odd-ordered  terms  in  the  equation 

.  .  .  .  r(v  +  77i  +  4)  .  . 

+ })(_„.).  (190) 

The  smallest  term  in  the  descending  series  occurs  where 

( m  -  I)2  -  vz 

—^ki - 1  (191) 

The  magnitude  of  the  smallest  term  is  given  by  the  Stirling  approximation 


_2Wm+j^  |\ve !- 

77771/  \m  -  V  —  =/ 


which  remains  bounded  as  v  -»  °°  if  the  ratio  between  m  and  v  exceeds  a  minimum 
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limit.  The  minimum  ratio  /3  satisfies  the  transcendental  equation 


whose  solution  is 


j8  =  1.543404638418208 

The  order  and  the  argument  are  bounded  by  the  limiting  relations 
m  =  pv  \z\  =  /iv 

for  which  the  ratio  g,  is  given  by  the  equation 

(Bs  -  1) 

a  =  — - -  =  0.4477432046943028 

^  2/S 


(193) 


(194) 

(195) 


(196) 


This  value  is  far  below  the  range  of  ratios  for  actual  use  of  the  asymptotic  series, 
which  can  be  cycled  until  there  is  no  change  in  the  sum. 

Programming 

SUBROUTINE  CBSSLJ  (AZ,  CN,  FJ) 

FORTRAN  SUBROUTINE  FOR  ORDINARY  BESSEL  FUNCTION  OF  FIRST  KIND 


The  real  and  imaginary  parts  of  the  argument  z  are  given  in  array  AZ,  and  the  real 
and  imaginary  parts  of  the  order  v  are  given  in  array  CN.  The  complex  Bessel  function 
is  computed  by  series  expansions  and  recurrence  relations.  Calls  are  made  to  Subroutine 
CGAMMA.  The  real  and  imaginary  parts  of  the  function  Jv(z)  are  stored  in  array  FJ. 
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MODIFIED  BESSEL  FUNCTION 

Analysis 

In  the  complex  plane  a  majority  of  the  Bessel  functions  are  multiple  valued.  A 
variety  of  Bessel  functions  for  real  argument  can  be  expressed  in  terms  of  the  modified 
Bessel  function  Kv(z)  for  complex  argument.  According  to  Equation  (8),  page  78  of 
Reference  24,  the  Hankel  function  of  the  first  kind  is  given  by  the  equation 

HvW(z)  =  -  -ie^vlriKv(-iz)  (197) 

7 r 

According  to  Equation  (5),  page  75  of  Reference  24,  the  Hankel  functions  are  related 
by  the  equation 

Hvw(zewi)  =  -  e~V7riHjs)(z)  (198) 

and  the  Hankel  function  of  the  second  kind  is  given  by  the  equation 

Hv(z\z)  =  -ie^vviKv(iz)  (199) 

7T 

Since  the  Bessel  functions  are  related  to  the  Hankel  functions  by  the  equation 

Jv(z)  =  ^HvM(z)  +  H™(z)  J  (200) 

the  Bessel  functions  are  given  by  the  equation 

Ju(z)  =  ~{esVmKv(iz)  -  e~zV7riKv{-iz)\  (201) 

and  since  the  Weber  functions  are  related  to  the  Hankel  functions  by  the  equation 

Yv{z)  =  ^HllM{z)-H^\z)])  (202) 

the  Weber  functions  are  given  by  the  equation 

Yv{z)  =  -  -{e^niKv(iz)  +  e~^vviKv{-iz)\  (203) 

7T  \  ) 

According  to  Equation  (18),  page  80  of  Reference  24,  the  modified  Bessel  function  /„(z) 
is  given  by  the  equation 

/„(*)  =  -kxze™)  -  e~v’niKv(z)}  (204) 

7 T  {  ) 

or  by  the  equation 

/„(z)  =  -  -\Kvize-™)  -  ev^Kv{z)\  (205) 

7T  (  ) 
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Elimination  of  7y(z)  from  these  two  equations  gives  the  equation 

Kv(zem)  +  Kv(ze~n%)  =  2  cos  vn  Kv(z)  (206) 

which  is  useful  for  extending  the  range  of  argument  beyond  the  usual  limits. 

The  Airy  function  is  defined  by  the  equation 

1  r 

Ai( x)  =  —  cos  ( xt  +  |f3)  dt  (207) 

^  Jo 

The  function  Ai(x)  is  expressed  by  the  equation 

Ai(x)  =  -(£*)®tfi(§**)  (208) 

tt  3 

and  its  derivative  Ai'(x)  is  given  by  the  equation 

Ai'(x)  = - i=  Kz{\ xz)  (209) 

7TV  3  3 

The  Airy  function  is  useful  for  the  approximation  of  Bessel  functions  when  the  argument 
and  the  order  are  large  and  nearly  equal.  They  have  applications  in  problems  of 
stationary  phase. 

On  the  positive  real  axis  Kv(z)  is  real  and  has  a  critical  region  where  it  cannot  be 
evaluated  easily  by  series  expansions.  It  is  a  natural  choice  for  a  rational  polynomial 
approximation  in  the  complex  plane.  The  argument  z  is  restricted  to  that  part  of  the 
complex  plane  from  which  the  negative  real  axis  is  excluded.  The  phase  of  z  is  limited 
to  the  range  -7T  to  +tt. 

If  the  argument  x  +  iy  satisfies  the  inequality 

xz  +  yz^\  (210) 


or  both  of  the  inequalities 


xz  +  yz  <  289  +  x  +  0.096t/2  g  0  (211) 

then  the  Bessel  functions  are  computed  from  the  ascending  power  series.  The  function 
Jy(z)  is  given  for  any  order  v  by  the  equation 


7„(z)  =  E 


1  „W+2m 


a=o  m!  r(v  +  m  +  1) 

whence  the  function  Kv(z)  is  given  for  fractional  order  by  the  equation 

t:  I_v(z)  -4(z) 

2  sin  vix 


(212) 


(213) 
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The  function  Kv(z)  is  given  for  integral  order  by  the  equation 


Kn(z)  =  (-  l)n+1  E 

m=0 


(I z)n+Zm 


m\(n  +  m)! 


7  +  log  (| z) 


1  m  i 

-  E  - 

2  “i  k 


■i  n+m 

-  E 

2  *=i 


*  n- 1 

+  w  E  (-D 


(n  —  m  — 
m! 


1)! 


(l2)-n+2m 


(214) 


where  7  is  Euler’s  constant.  The  summations  are  continued  until  there  is  no  change 
in  the  associated  value  of  7„(z). 

If  the  argument  x  +  iy  satisfies  both  of  the  inequalities 

1  <  xz  +  yz  <  289  +  x  +  0.096  y2  >  0  (215) 


then  the  Bessel  function  is  computed  from  the  rational  approximation.  When  the 
converging  factor  is  approximated  by  poles  and  residues,  then  the  Bessel  function  is 
expressed  by  the  equation 


Kv{z) 


1  +  E 


Zi  *  -  <5* 


(216) 


for  which  the  positions  6i  and  the  residues  ei,  for  orders  0,  |,  §,  1,  are  listed  in  Tables 
VII  through  X. 

If  the  argument  x  +  iy  satisfies  the  inequality 


xz  +  yz>2&9 


(217) 


then  the  Bessel  functions  are  given  by  descending  power  series.  According  to 
Equation  (4),  page  172  of  Reference  24,  the  integral  representation  of  the  function 
Kv(z)  is  given  by  the  equation 

(2ie) 

which  may  be  transformed  through  the  substitution 

t  =  1  +  -  (219) 

z 


to  give  the  equation 


Kv{z) 


7T  \2 


/IOC 

T) 

2/  VO 


e-V'-tf  1  +  U 


,2z)  T(v  +  2 

Expansion  of  the  integrand  with  the  binomial  series 


V 


n^P-E1  r(y  +  i) 

2 z)  m=0  m\Y{v  -  m  +  \)\2z  1 


i>  +  |) 


_ (±\* 

(N-  l)!r(i/-JV  +  |)\2z/ 


2z) 


du 


(220) 


\yt) 


V-AT-I 


(1  -  f)*-1  dt 


(221) 
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and  term  by  term  integration  leads  to  the  equation 


1 

Kv^  U*)  6  ™,\r(v  -  m  +  |)(2z)r 


^  T{v  +  m  +  i) 


(AT  -  l)!(2z) 


2^/ 


(1  -  t)N~1dt  du 


(222) 


in  which  the  series  is  asymptotic.  The  function  is  given  by  the  equivalent  equation 

(223) 


KM  '  (§)V 


tf-1  \ 

E  ^m(z)  +  t/Ar(z)AA,(f7) 

m=0  / 


where  the  asymptotic  terms  Um(z)  are  defined  by  the  equation 


vM  = 


(Vz  _  (i)*)(y*  _  (|)8)-(^  -  (m  -  |)a)(^  -  (m  -  |)2) 
m!(2z)m 


(224) 


and  Ay(<j)  is  a  converging  factor.  According  to  the  formulae  on  page  547  of  Reference  1, 
the  converging  factor  hN{o)  is  given  in  terms  of  the  deviation 

a  =  \/xz  +  yz  -  gAT  (225) 

and  the  tangent 


tan  6  = 


y 

x  +  Va;2  +  yz 


(226) 


by  the  equation 


A*/  — 


a  1  \  sec20 


+  ~  +  8 


+  l 


N 


a  a  3 

—  +  —  + - 

2  8  32 


Bvz  1 3a  5  ,  , 

+  —  +  —  tan20 
8  32 


sec20 

BN 


a  a  l 
—  -f  —  +  — 
2  4  32 


32 


tan20 


sec2# 

Nz 


sec  z6 
Nz 

tan  6 


(227)  , 


With  this  converging  factor  and  double  precision  it  is  possible  to  obtain  better  than 
single  precision  values  for  use  in  rational  approximation. 

Otherwise,  summation  of  the  asymptotic  series  is  continued  until  there  is  no  change 
in  the  sum  or  until  m  =  35. 

After  the  Bessel  functions  of  lowest  order  have  been  computed,  the  Bessel  functions 
of  higher  order  can  be  generated  with  the  aid  of  the  recurrence  equation 


Kv+1(z)-Kv_1(z)  =  —Kv(z) 
z 


(228) 


A  rounding  error  eM  is  introduced  during  the  /zth  cycle  of  recurrence  and  persists 
through  the  subsequent  cycles.  The  error  can  be  expressed  as  a  linear  combination 
of  the  functions  (—l)vIv(z)  and  Kv(z)  both  of  which  obey  the  same  recurrence  equations. 
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In  view  of  the  recurrence  equation 

lv(z)Kv+i(z)  +  Iv+1(z)Kv(z)  =  -  (229) 

z 

the  persisting  error  in  the  vth  cycle  is  given  by  the  expression 

z  j(-l)^_1(2)/u(2)+V1(2)^(2)JeM  (230) 

Insofar  as  the  function  Kv(z )  increases  with  order  faster  than  Iv(z),  the  relative 
rounding  error  remains  bounded. 

Programming 

SUBROUTINE  CBSSLK  (AZ,  CN,  FK) 

********************************************************************************************* 
FORTRAN  SUBROUTINE  FOR  MODIFIED  BESSEL  FUNCTION  OF  SECOND  KIND 
********************************************************************************************* 

The  real  and  imaginary  parts  of  the  argument  2  are  given  in  array  AZ,  and  the 
order  v  is  given  in  CN.  The  complex  Bessel  function  of  one— third  integral  order  is 
computed  by  series  expansions,  rational  approximations,  and  recurrence  relations. 
The  real  and  imaginary  parts  of  the  function  Kv(z)  are  stored  in  array  FK. 
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DISCUSSION 


The  subroutines  for  this  report  were  programmed  originally  for  the  Naval  Ordnance 
Research  Calculator.  When  the  NORC  was  about  to  be  demolished,  the  subroutines  were 
converted  to  a  hybrid  FORTRAN  which  would  run  with  the  STRETCH  compiler.  Register 
storage  in  the  NORC  had  the  same  function  as  the  accumulator  in  the  STRETCH,  but 
there  is  no  way  to  refer  to  the  accumulator  through  FORTRAN,  so  the  register  storage 
was  simulated  by  a  storage  location  in  memory.  This  contributed  to  a  two— fold  loss 
of  efficiency  which  has  been  observed  with  the  hybrid  FORTRAN-  The  subroutines  have 
been  reprogrammed  from  their  formulations  directly  into  FORTRAN  in  order  to  improve 
efficiency  and  accuracy. 

In  the  present  work  the  approximation  of  converging  factors  was  accomplished  in 
two  steps,  of  which  the  first  was  an  adjustment  of  the  rational  approximation  and  the 
second  was  a  conversion  to  a  polar  expansion. 

Future  development  of  expansions  in  poles  and  residues  probably  will  take  the  form 
of  transformations  to  Stieltjes  integrals  which  then  can  be  approximated  by  Laguerre 
integration.  Given  any  asymptotic  expansion,  the  integrand  of  the  Stieltjes  integral 
can  be  so  constructed,  from  a  series  of  Laguerre  functions,  as  to  have  the  same 
asymptotic  expansion,  whence  Gauss-Laguerre  quadrature  would  lead  to  an 
approximation  in  poles  and  residues. 

Another  approach  to  the  approximation  of  functions  follows  from  the  fact  that  the 
real  and  the  imaginary  parts  of  an  analytic  function  are  solutions  of  Laplace’s  equation. 
A  sourcewise  representation  of  the  solution  of  Laplace’s  equation  consists  of  poles 
with  residues. 


CONCLUSION 

A  number  of  functions  can  be  approximated  successfully  by  sets  of  poles  with 
residues.  The  functions  first  were  approximated  as  rational  functions,  then  were 
converted  to  the  equivalent  singular  functions. 
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APPENDIX  A 


ALGORITHMS 


METHOD  OF  APPROXIMATION 


Rational  Functions 

A  double  precision  complex  routine  was  available  on  NORC  for  the  rational 
approximation  of  a  function  within  a  half  circle  of  unit  radius.  Input  consisted  of 
values  of  the  argument  2  and  values  of  a  function  R(z)  so  closely  spaced  that 
interpolation  gave  full  accuracy.  The  program  constructed  approximations 


R(z) 


?>(z) 

q(z) 


(1) 


of  progressively  higher  order  until  a  tolerance  was  met.  Output  consisted  of  the  values 
of  the  coefficients  of  the  polynomials  p(z)  and  q{z). 

In  the  rational  function  of  the  nth  degree  there  were  2n  +  1  complex  coefficients, 
which  were  specified  by  the  2n  +  1  complex  values  of  the  function  to  be  approximated, 
at  2n  +  1  nodes  on  the  complex  contour.  The  origin  at  2  =  0  and  the  corners  at  2  =  ±i 
were  included  among  the  nodes,  but  since  there  must  be  an  odd  number  of  nodes, 
the  point  at  2  =  ±1  could  not  be  a  node. 

The  function  R(z)  was  specified  to  be  symmetric  with  respect  to  the  real  axis  and 
the  polynomials  p(z)  and  q(z)  had  real  coefficients.  Values  of  the  function  on  the 
positive  quarter  circle  were  sufficient  to  establish  the  real  coefficients.  The  method 
of  inverted  differences  was  used  to  determine  the  coefficients.  The  relative  error  £(2) 
was  computed  for  each  datum  of  input. 

The  coordinates  of  the  nodes  were  adjusted  to  make  the  maxima  of  error  between 
nodes  more  uniform.  The  coordinate  of  the  fcth  node  on  the  imaginary  axis  was  the 
distance  yk  from  the  origin  and  the  coordinate  of  the  fcth  node  on  the  circular  arc 
was  the  angle  6k  with  respect  to  the  real  axis. 

Inasmuch  as  the  error  was  analytic  within  the  complex  contour  it  could  be  factored 
into  a  polynomial  which  had  roots  at  the  nodal  points  on  the  contour.  The  error  £(2) 
thus  could  be  expressed  by  the  equation 

k=2N+l 

e{z)=A{z)  []  (z-z*)  (2) 


where  A(z)  is  an  analytic  function  and  zk  is  the  fcth  nodal  point.  Insofar  as  A{z)  was 
nearly  constant,  an  increment  in  error  Ae  was  expressed  in  terms  of  displacements 
of  nodes  hzk  by  the  equation 


*=1  (z-z*) 

The  increment  A|e|  in  the  amplitude  of  error  |e|  was  given  by  the  equation 

Azk 


(3) 


A|£ 


K  =  cLN+ 1  - 

=  -  m  £  ; 

*= 1  L ' 


‘■k  +  A2£ 


2*  Z*  -  Z? 


(4) 
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On  the  imaginary  axis  the  increment  A zk  was  given  in  terms  of  distance  A yk  by  the 
equation 

A  zk  =  iAyk  (5) 

and  on  the  unit  circle  the  increment  A zk  was  given  in  terms  of  angle  A 8k  by  the 
equation 


A  zk  =  iel6,c  A6k 


(6) 


Equalization  of  two  successive  maxima  of  error  was  expressed  by  the  equation 


(?) 


where  |e|m  was  the  amplitude  of  error  at  the  mth  maximum  of  error.  Substitutions 
from  Equations  (4),  (5),  (6)  in  Equation  (7)  led  to  a  system  of  simultaneous  equations 
which  could  be  solved  by  a  standard  method  to  give  increments  A yk,  A 8k. 

After  the  maxima  of  error  were  uniform  along  the  imaginary  axis  and  were  uniform 
along  the  circular  arc,  the  number  of  nodes  on  either  line  was  increased  according 
to  whichever  line  had  the  larger  maxima  of  error.  The  coordinates  on  the  line  were 
interlineated  to  make  space  for  the  added  node. 

This  algorithm  for  approximation  in  a  finite  zone  within  a  half  circle  of  unit  radius 
was  adapted  to  the  problem  of  approximation  in  the  infinite  zone  outside  a  circle  of 
unit  radius  when  the  variable  z  was  replaced  by  the  reciprocal  z  in  the  choice  of 
argument. 

Equivalent  Poles 

A  ratio  between  two  power  polynomials  can  be  expressed  as  the  sum  of  singular 
terms  if  the  degree  of  the  numerator  is  less  than  the  degree  of  the  denominator.  The 
singular  terms  represent  poles  which  are  situated  at  the  roots  of  the  denominator. 
Let  the  ratio  R(z)  be  expressed  by  the  equation 


sw.ja.sizivjfr: 

q{z)  c0  +  •••  +  cnzn 


(8) 


where  cn  is  equal  to  unity.  Let  the  roots  of  q(z)  be  z1 . zn.  It  is  possible  to  take  out 

a  term  which  contains  only  zn,  as  expressed  by  the  equation 


R(z) 


+  a'p  +  -  +  a’n-zz71  8 

z-zn  c’o  +  -  + 


0) 


where  rn  is  a  constant,  a.g,  ...,  a'n_z  are  the  coefficients  of  a  polynomial  of  lower  degree, 
and  Co . c^j  are  the  coefficients  of  the  quotient 


Cq  + 


+  c' 


cn  +  •••  +  c„zT 


(10) 


In  order  to  evaluate  the  coefficients,  let  Equation  (9)  be  reduced  to  a  common 
denominator,  whence  the  numerator  becomes  p(z)  and  the  denominator  becomes  q(z). 
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When  the  terms  of  the  same  degree  in  z  are  equated  on  both  sides,  then  the  result 
is  the  following  system  of  equations. 


an  =  cnr„ 


a,  =  c,r„ 


a'hZn  +  a'k-l  > 


(11) 


“n-1  =  c'n-irn  +  a'n_z  J 

If  these  equations  are  multiplied  in  succession  by  the  factors  1 . znk . znn_1  and 

are  added,  then  the  last  two  terms  on  the  right  cancel  in  pairs,  and  the  terms  which 
remain  lead  to  the  equation 


an  +  —  +  a__,  z_ 


rn  = 


n  I  o'  Tl-1 

Thus  it  is  possible  to  evaluate  rn  by  a  straightforward  computation. 
Rearrangement  of  the  Equations  (11)  leads  to  the  chain  of  equations 


^n-2  ^n— 1  ^n-l^n 


(12) 


a'k-i  =  a*  ~  c'krn  +  a'kZr 


(13) 


ao  =  at  -  c\rn  +  a\zn 

which  can  be  solved  successively  to  obtain  the  coefficients  a'n_z . a'0  of  the  reduced 

polynomial.  The  ratio  R{z)  thus  can  be  expanded  recursively  into  the  sum 


R(z) 


ri  rn 

— —  +  + - - 


z  -  z. 


z  —  z_ 


(14) 


where  zi,  ....  zn  are  the  positions  of  poles  and  r1,...,rn  are  the  residues  of  the  poles. 
In  the  application  of  this  algorithm  to  rational  approximations,  the  ratio 


p(z)  a0  +  •••  +  anzn  an  +  ■■■  +  a  azT 


(15) 


q(z)  c0  +  ---  +  cnzn  cn  +  ■■■  +  c0zn 
can  be  resolved  into  the  sum  of  the  two  terms 

1  +  (gn  -  £n)  +  •"  +  K  -  cjz71"1 
Crv  +  •••  +  C0Zn 

because  a0  =  c0,  and  then  the  second  term  can  be  expanded  into  singular  terms.  The 


(16) 


3 


factor  z  can  be  included  in  the  expansion  to  give  the  equation 


p(z) 

q(z) 


m 


which  contains  the  pole  at  the  origin,  while  the  residues  are  those  for  the  expansion 
of  the  second  term  in  Expression  (16). 
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APPENDIX  B 


RATIONAL  APPROXIMATIONS 


TABLE  I 


COMPLEX  EXPONENTIAL  INTEGRAL 

Ei(z) 

POLES 

Positions  Residues 


0.000000000000000 

E 

00 

8.501565161210931 

E— 03 

3.111059570865283 

E- 

02 

5.050374658490585 

E— 02 

1.036612605391116 

E- 

■01 

8.368173689564071 

E-02 

2.165323352445536 

E- 

■01 

1.070475824176067 

E— 01 

3.699314279601916 

E- 

-01 

1.204247190294619 

E— 01 

5.667662599905892 

E- 

-01 

1.250966315822293 

E— 01 

8.140420663247483 

E- 

-01 

1.223144352246853 

E-01 

1.123842475408128 

E 

00 

1.126214175539071 

E— 01 

1.514004781485123 

E 

00 

9.634194073925819 

E-02 

2.008867950322836 

E 

00 

7.473984227575106 

E-02 

2.640524118235915 

E 

00 

5.085961359534413 

E-02 

3.450984499333923 

E 

00 

2.908227067736279 

E-02 

4.495833607632020 

E 

00 

1.322016405301009 

E-02 

5.850582634098224 

E 

00 

4.438029398290675 

E-03 

7.622735014633804 

E 

00 

9.926124789875764 

E— 04 

9.978145015845783 

E 

00 

1.265797951120111 

E— 04 

1.321220648964083 

E 

01 

7.021509082533505 

E-06 

1.803229483760214 

E 

01 

9.102815325646319 

E-08 

TABLE  II 


COMPLEX  FRESNEL  INTEGRAL 

E(z) 

POLES 


Positions 


0.000000000000000 

E 

00 

2.086058560134765 

E- 

■02 

8.298069404956873 

E- 

-02 

1.854216533260787 

E- 

-01 

3.279634793823607 

E- 

-01 

5.126752799128284 

E- 

-01 

7.454129580451047 

E- 

-01 

1.036950674182965 

E 

00 

1.403780612554370 

E 

00 

1.868916622140010 

E 

00 

2.463148305239293 

E 

00 

3.227193837373523 

E 

00 

4.215343482800130 

E 

00 

5.501788731515490 

E 

00 

7.192589666831019 

E 

00 

9.451702080764080 

E 

00 

1.257107183147839 

E 

01 

1.724835372163339 

E 

01 

Residues 

8.157230833240962  E-02 
1.592852852534368  E-01 
1.485816256144991  E-01 
1.332196708362453  E-01 
1.156903928789572  E-01 
9.785809594475354  E-02 
8.059088342976243  E-02 
6.402045386098722  E-02 
4.814452427678847  E-02 
3.335406584732945  E-02 
2.055480994701934  E-02 
1.078474038875057  E-02 
4.556348922142192  E-03 
1.439844581389254  E-03 
3.070561398341705  E-04 
3.781565411685414  E-05 
2.051735096161211  E-06 
2.635648236827474  E-08 


TABLE  III 


DOUBLE  EXPONENTIAL  INTEGRAL 

Di(z) 

POLES 


Positions 

Residues 

0.000000000000000 

E 

00 

-3.469117335358920 

E— 04 

4.195566783742928 

E- 

-02 

-6.037877324617450 

E— 03 

1.175336616486651 

E- 

-01 

-1.524613059492494 

E-02 

2.285602374559868 

E- 

-01 

-2.105821698272908 

E-02 

3.756673501612400 

E- 

-01 

-1.718942087207540 

E-02 

7.915948462766724 

E- 

-01 

3.143234670330322 

E-02 

1.075468896230582 

E 

00 

7.508985315669721 

E-02 

1.426592080308411 

E 

00 

1.246897878072599 

E-01 

1.862905549523766 

E 

00 

1.685790750900348 

E— 01 

2.407300095098562 

E 

00 

1.917150806995113 

E-01 

3.088540356075242 

E 

00 

1.826007945508359 

E-01 

3.942776051552587 

E 

00 

1.423456743071465 

E-01 

5.015931965439807 

E 

00 

8.748622224193267 

E-02 

6.367591807486506 

E 

00 

4.021750832884253 

E-02 

8.077761930960550 

E 

00 

1.285750056801799 

E-02 

1.025989611388872 

E+01 

2.576737825984412 

E— 03 

1.308967684226102 

E+01 

2.759550037843490 

E-04 

1.688321690859160 

E+01 

1.191393155171215 

E— 05 

2.240832409417126 

E+01 

1.072929801993859 

E— 07 

TABLE  IV 


DOUBLE  FRESNEL  INTEGRAL 

D(z) 

POLES 


Positions  Residues 


0.000000000000000 

E 

00 

1.946702176629290 

E-03 

6.158254330058267 

E- 

•03 

-5.602038332003368 

E-03 

5.577129558799173 

r_ 

•02 

-7.734209998293494 

E-03 

2.613294506565088 

E- 

-01 

1.923681980789794 

E-02 

4.197557282462471 

E- 

-01 

4.685483890266400 

E-02 

6.191129560306540 

E- 

01 

7.868054439865047 

E-02 

8.663690369835087 

E- 

■01 

1.103380471995033 

E-01 

1.172669313738075 

E 

00 

1.367629763844977 

E-01 

1.553763009533203 

E 

00 

1.519362601026446 

E-01 

2.030480287272866 

E 

00 

1.502574885873442 

E-01 

2.629550653406246 

E 

00 

1.296972574495706 

E-01 

3.385002037633292 

r 

00 

9.482725744458265 

E-02 

4.340433251697753 

E 

00 

5.639680002398576 

E-02 

5.552727526730474 

E 

00 

2.585885948785604 

E-02 

7.098490802164333 

E 

00 

8.504554154209225 

E-03 

9.086373817410945 

E 

00 

1.812509189092737 

E-03 

1.168418085479360 

E+01 

2.144626967496620 

E— 04 

1.519234066877839 

E+01 

1.074839393817003 

E— 05 

2.033584655600194 

E+01 

1.219304804726968 

E— 07 

TABLE  V 


EXPONENTIAL  EXPONENTIAL  INTEGRAL 

k{z) 

POLES 

Positions  Residues 


0.000000000000000 

E 

00 

2.372861283136832 

E- 

-02 

8.541132106687596 

E- 

-02 

1.852766272820592 

E- 

-01 

3.237415266166881 

E- 

-01 

5.030454603812668 

E- 

-01 

7.288066075871879 

E- 

-01 

1.011227701028717 

E 

00 

1.365984481712494 

E 

00 

1.815101390389295 

E 

00 

2.388247019554190 

E 

00 

3.124905320088117 

E 

00 

4.078024898944454 

E 

00 

5.320335455548645 

E 

00 

6.956243072905789 

E 

00 

9.147599025470312 

E 

00 

1.218290743885443 

E 

01 

1.675113119698726 

E 

01 

3.495172589268265  E-02 
1.358491059258968  E-01 
1.588505815522956  E-01 
1.530014345354354  E-01 
1.345207528564613  E-01 
1.119130516196707  E-01 
8.920083866561896  E-02 
6. 792272054720666  E-02 
4.867231978872109  E-02 
3.201709765322659  E-02 
1.870089650211108  E-02 
9.297084144278655  E-03 
3.726047631610865  E-03 
1.119895375598227  E-03 
2.280574968723531  E-04 
2.695962277814528  E-05 
1.412554302243009  E-06 
1.763523268088063  E-08 


TABLE  VI 


EXPONENTIAL  FRESNEL  INTEGRAL 
m(z) 


POLES 


Positions 

0.000000000000000  E  00 
1.882000761332742  E-02 
7.650297830721689  E-02 
1.731094073268224  E-01 
3.087247749376508  E-01 
4.852407450425483  E-01 
7.079315023126571  E-01 
9.866948545085860  E-01 
1.336964217293807  E  00 
1.780586961924434  E  00 
2.347109262385231  E  00 
3.075935359899532  E  00 
4.019969426341625  E  00 
5.251927395269674  E  00 
6.876223046863616  E  00 
9.054707180068234  E  00 
1.207577867202912  E  01 
1.662778013418871  E  01 


Residues 

1.171632185997974  E-01 
2.084876153495584  E-01 
1.720712026818202  E-01 
1.374456548793750  E-01 
1.071559279197100  E-01 
8.191163553714233  E-02 
6.130560309120205  E-02 
4.448272993284047  E-02 
3.070308658007574  E-02 
1.961788115885990  E-02 
1.120676960634960  E-02 
5.480694573733344  E-03 
2.171717980701646  E-03 
6.482424452554320  E-04 
1.316220788022128  E-04 
1.556852300136330  E-05 
8.187688247021873  E-07 
1.029295020822967  E-08 


TABLE  VII 


BESSEL  FUNCTION  OF  ORDER  ZERO 
K0(z) 

POLES 

Positions 


0.000000000000000  E  00 
•1.648995051422117  E-02 
■7.186218800685365  E-02 
1.670868781248656  E-01 
■3.025822502194688  E-01 
■4.806139452459267  E-01 
■7.070752393578979 1  E-01 
9.929957905395160  E-01 
1.355839256125922  E  00 
1.821059078991320  E  00 
2.424821753108787  E  00 
3.219566557087496  E  00 
4.286580772483836  E  00 
5.770228167981279  E  00 
8.013712609525260  E  00 


Residues 

0.000000000000000  E  00 
-4.809423363874473  E-03 
-1.313662003477595  E-02 
-1.948438340084579  E-02 
-2. 1994B9000320033  E-02 
-2.093966256765194  E-02 
-1.746002684586503  E-02 
-1.279378133620848  E-02 
-8.052344217965918  E-03 
-4.158173750027601  E-03 
-1.643177387479224  E-03 
-4.491755853147087  E-04 
-7.285947655740069  E-05 
-5.382652306582855  E-06 
-9.937790480362892  E-08 


TABLE  VIII 


BESSEL  FUNCTION  OF  ORDER  ONE-THIRD 


Ki(z) 


POLES 


Positions 

0.000000000000000  E  00 
1.436346451938738  E-02 
6.829706137621937  E-02 
1.626145943132923  E-01 
2.974538543947166  E-01 
■4.749052692100660  E-01 
7.007380209350545  E-01 
9.858961485754878  E-01 
1.347791024780270  E  00 
1.811850277214040  E  00 
2.414232693109992  E  00 
3.207380096917656  E  00 
4.272584301588853  E  00 
5.754199732132192  E  00 
7.995346895730257  E  00 


Residues 

0.000000000000000  E  00 
-3.308130422855452  E-03 
-7.992541489780344  E-03 
-1.114953951149366  E-02 
-1.214944767056973  E-02 
-1.131050261572697  E-02 
-9.290616383535283  E-03 
-6.738338829864222  E-03 
-4.212029053420502  E-03 
-2.165634872117979  E-03 
-8.537655388612388  E-04 
-2.332000490440423  E-04 
-3.784577154525431  E-05 
-2.800396179763318  E-06 
-5.183943839708712  E-08 


TABLE  IX 


BESSEL 


Positions 

0.000000000000000 

-1.002408161647787 

-5.997870547395276 

-1.516343542993466 

-2.845368398965238 

-4.603218236550183 

-6.844133505234451 

-9.675114349534103 

-1.326876006320234 

-1.787857730177285 

-2.386589640411052 

-3.175517091435769 

-4.235942093386486 

-5.712193661037926 


FUNCTION  OF  ORDER  TWO-THIRDS 
Kz{z) 

3 


POLES 

Residues 


E  00 

0.000000000000000 

E  00 

E— 02 

8.356269613036070 

E-03 

E— 02 

1.421589792080583 

E-02 

E— 01 

1.661492960434507 

E-02 

E— 01 

1.632135891550809 

E-02 

E-01 

1.421532926606566 

E-02 

E— 01 

1.116165182456759 

E-02 

E-01 

7.846521235416155 

E-03 

E  00 

4.800776505229076 

E-03 

E  00 

2.433841308314365 

E-03 

E  00 

9.515251496981341 

E-04 

E  00 

2.589249109738377 

E-04 

E  00 

4.201905951784872 

E-05 

E  00 

3.118825106795326 

E-06 

-7.947172476979304  E  00 


5.808362400385792  E-08 


TABLE  X 


BESSEL  FUNCTION  OF  ORDER  ONE 
K,{z) 

POLES 

Positions 


0.000000000000000  E  00 
5.577424298795054  E-03 
4.991129441724757  E-02 
1.374409116523968  E-01 
2.672337847105657  E-01 
4.403801668086817  E-01 
6.618136148725406  E-01 
9.418610776650166  E-01 
1.297541304683261  E  00 
1.754076967198164  E  00 
2.347552998822763  E  00 
3.130413326891964  E  00 
4.183971205637291  E  00 
5.652517992149936  E  00 
7.878639598106769  E  00 


Residues 

0.000000000000000  E  00 
7.538057792005914  E-02 
7.122935374034643  E-02 
6.331162242281997  E-02 
5.282402645233011  E-02 
4.133053594414916  E-02 
3.013505739475096  E-02 
2.010434395927201  E-02 
1.185522230680744  E-02 
5.860555109560099  E-03 
2.254651482673253  E-03 
6.081730415363355  E-04 
9.842155506257467  E-05 
7.321390930380890  E-06 
1.372796673846658  E-07 


APPENDIX  C 


ERRORS  OF  APPROXIMATION 


FUNCTION 

Ei(z) 

E(z) 

Di(z) 

D(z) 

k(z) 

m(z) 

K0  (z) 

Ki  {z) 

3 

Kz{z) 

3 

K,{z) 


MAXIMUM  RELATIVE  ERROR 

ERROR 

Axis 

1 1.6x  10~16 

4.1xl0“le 

5.1xl0“16 

4.8xl0“16 

2.8x  10~16 

1.7xl0“ls 

3. 1  x  10-16 

1.7xl0“16 

2.0xl0-16 

5.1X10"16 


Arc 


4.4X10-16 


1.8xl0-le 


2.2x  10~16 


2.0xl0"16 


1 ,4x  10~16 


0.9x  10_le 


1.3xl0"16 


0.7x  10~16 


0.8x  10~16 


2.1xl0"16 
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Internal 

D 

K 

KR 

KA 

KC 

KE 

KG 

KO 

KP 

KW 

KXD 

KXF 

KXH 

KXS 

KXU 

KXW 

MAL 

File 


1 

1 

1 

1 

1 

1 

1 

1 

3 

1 

2 

1 
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1 

2 

1 

5 
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